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I. 



INTRODUCTION 



A. PURPOSE 

Since the end of the Cold War, there has been a paradigm shift in detecting, 
tracking, classifying and ultimately neutralizing subsurface threats. In its most general 
form, the US Navy now sees a more probable warfare scenario in coastal waters as 
opposed to the Blue Water confrontations occurring earlier last century. The major threat 
of the once large Soviet Nuclear Submarine Force has now turned into a small-country 
coastal diesel submarine threat. Indeed, the very name for this warfare area changing 
from Anti-submarine Warfare to Undersea Warfare implies a change in our subsurface 
warfare strategy. 

With the shift in subsurface threats come different challenges in detection and 
classification. No longer can acoustic detection systems exploit deep sound channels on 
a relatively noisy nuclear submarine threat. In this new era where the threat of the once 
large nuclear fleet from the Former Soviet Union is limited at best, a new foe has 
emerged in coastal waters. Not only do Naval USW systems have to contend with 
increased challenges such as bottom bounce and added noise due to a higher density of 
neutral shipping, but the platforms of choice are the quiet diesel submarines while on 
battery power. In addition to this diesel threat in a Littoral Warfare environment, mines 
have become more of a threat. Because of their high destruction-to-cost ratio, laying 
mines to thwart US global interests is a cost-effective way to inflict maximum damage 
with little cost. Indeed, the US Navy has suffered more battle damage to their ships due 
to mines than any other subsurface threat since the Korean War. Mines are cheap, easy to 
deploy and hard to detect. To many countries, mines are a viable substitute for an 
otherwise weak coastal defense. 

The nature of Littoral Warfare possesses unique challenges to detecting and 
classifying subsurface threats. Since the Walker Treason in the middle 19S0's, enemy 
acoustic signatures have been suppressed to an all-time low. One way to regain a 
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considerable edge in Undersea Warfare is to investigate more robust methods of detecting 
enemy acoustic signatures within heavy background noise. 



B. GOALS 

This thesis is part of a larger project that addresses the issue of improving the 
tools available to USW systems operators of surface and subsurface platforms. 
Currently, USW systems operators use variations of a spectrogram as their main visual 
display to monitor passive acoustic signatures. This type of analysis is good for signals 
that are long in duration. Existing displays show a long-term history of acoustic data, 
which is good for determining underwater activity over a long period of time. Using the 
same techniques, short-duration signals (or transients), however, may go unnoticed. For 
example, if an operator is either fatigued or focusing on long-term trends, he might 
overlook a very slight indication that a transient signal a fraction of a second in duration 
had just occurred. 

The thesis goal is to create more robust visual and audible features at the 
operator’s disposal to enhance the detection of transient signatures. Specifically, the 
thesis will develop and investigate three signal processing procedures based on Wavelet 
and Multiresolution analyses, which will attempt to detect transient signals currently 
overlooked by common signal processing techniques. 

This thesis will also investigate an automation procedure, which will alarm the 
operator of such transient signals. 



C. METHODOLOGY 

It is the main thrust of this thesis to detect transient signals in a noisy acoustic 
environment. One way to investigate the improvement of detecting desired acoustic 
signatures imbedded within a noisy environment is by studying the statistics and nature 
of the signals. Although littoral warfare might be difficult in terms of background 
shipping and silent threats, one positive aspect is that acoustic data is easily accessible 






because the surveillance area is small and there exist well-known shipping lanes and 
underwater topographies. 

Current methods of transient detection use time-frequency analysis and stochastic 
modeling. Although these methods have been proven quite useful, they have their 
limitations. This thesis will investigate the root causes of these weaknesses and explore 
some ideas to circumvent them. Ideas such as Filter Banks using Wavelet Transforms 
and Multiresolution Filtering just might be the missing tools to enhance transient 
detection. 

Three methods have been developed to detect transient signals in noise. One of 
these methods is based on Wavelet Analysis and the other two based on an ARMA model 
and multiresolution processing. A standard data set is used to compare the new signal 
processing methods against well-established methods. From this data, conclusions will 
be drawn and recommendations will be made for further improvement of acoustic 
transient detection systems. 



D. BENEFITS 

The confidence of a decision is directly proportional to the amount of concurring 
evidence. Consider the method of detecting submarines using the classic method of 
Target Motion Analysis (TMA). Since TMA uses nothing but passive techniques, 
determining a target’s bearing, range, course and speed can be challenging. What 
strengthens this method of detecting underwater weapons systems is the redundancy of 
evidence. Information that does not show up on a Dead Reckoning Trace (DRT) plot 
might be found on the time-frequency table or on a Moboard plot. Given a well-behaved 
target, the success of TMA relies on the fact that dynamic changes will be detected in 
different domains. For example, a shift in frequency may be the first detection of a 
change in course or speed, but probably won’t be the first method to estimate a target’s 
range or course. The different techniques in TMA are independent in the sense that they 
do not rely on similar assumptions. These methods each provide a piece of a puzzle, 
which when put together, can accurately determine particular warfare scenario. 
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Like the set of tools used in TMA, an acoustic transient detection system needs to 
be developed which takes the same fact-layering approach in order to make an informed 
decision. This thesis does not attempt to improve on well-established signal processing 
techniques such as time-frequency analysis. Rather, the goal of this thesis is to create 
additional layers (based on fundamentally different assumptions) and add to the fact- 
finding set in order to enhance the robustness of detecting transient signals. 

The benefits of this thesis will be the addition of new tools, which will enhance a 
transient detection system. One of the new signal processing techniques may be the 
necessary tool needed to detect some transient signals. 

Another benefit to the undersea warfighting capability of a ship is to discuss 
methods to enhance the automation of detecting valid transient signals. 
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II. 



CURRENT TRANSIENT SIGNAL DETECTION TOOLS 



Although current methods of detection can successfully identify transient signals, 
there is a threshold floor or frequency bandwidth to which signals become 
indistinguishable from background noise. Of those cases where current detection 
schemes miss a signal, we will show that the new methods discussed in this thesis will 
have a strong indication. This chapter establishes strengths and weaknesses of existing 
signal detection methods both in the time and frequency domains. The intent is not to 
produce an all-inclusive list of current signal processing techniques, but rather to 
compare new methods against existing ones. In this way, we hope to show that these new 
techniques will not replace current ones, but rather support them as they add to the 
robustness of the whole Acoustic Transient Detection System. 



A. ESTABLISH A BENCHMARK SIGNAL 

In order to meet the goal of enhancing transient signal detection, a comparison 
between established methods and new methods presented in this thesis must be 
performed. To accomplish this evaluation, a single representative signal will be used as a 
benchmark. The model signal is the sound of a fan blowing for six seconds sampled at 
11.025 kHz. Starting after the third second, a series of 11 tapping sounds (each 
approximately .1 seconds in duration) will occur. These taps will be part of the analysis 
on transient resolution. The first three seconds will be used to study the effectiveness of 
all transient signal detection techniques by varying the amplitude, duration, and 
frequency content of an added transient signal. Figure II- 1 shows a block diagram of the 
filtering process of the signal to detect the transient noise. Figure II-2 is a pictorial 
representation of the benchmark comparison setup based on the filtering process 
described in Figure II- 1. 
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Figure II-l Block Diagram of Whitening Filter Process of Benchmark Signal 
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Figure II-2 Whitened and Filtered Acoustic Signal Benchmark (z[n]) 
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1 . 



What Kind of Noise? 



It is widely held that since acoustic ambient noise is comprised of many random 
sources, it can be assumed to be Gaussian [2]. This assumption is based on the Central 
Limit Theorem, which states that as the number of random variables in a sequence 
becomes large, the probability distribution of the sum of the sequence approaches a 
Gaussian random variable [4], In the case of shallow water, however, ambient ocean 
noise coming from the Infrasonic Band and the Low Sonic Band may dominate, and 
therefore may compromise the previous assumption. 

The Infrasonic Band comprises frequencies between 1 Hz and 20 Hz. This is the 
frequency band that contains blade rate fundamental frequencies of power-driven vessels. 
The Low Sonic Band has a frequency range between 20 Hz and 200 Hz, which is mostly 
caused by distant shipping, especially in the frequencies around 100Hz [2], 

At any rate, noise can be defined simply as any part a signal that you don’t want. 
For the specific case of detecting transient signals in shallow water, it is a good 
assumption that the above-mentioned ambient noise bands dominate the noise model. 

Take for example the scenario of attempting to detect a small ship in shallow 
water laying mines among a fleet of fishing vessels. It is desirable to detect a certain 
transient signal like a splash or machinery among many acoustic sources with similar 
frequency content. Since we know that mines are deployed in heavily traveled 
waterways, you can assume that the background noise (defined as the part of the signal 
you don’t want) will have a dominant low frequency content. The model therefore, will 
not be white noise, but colored noise since it has uneven frequency content. 

The colored noise assumption is reflected in the Benchmark Signal. Much like a 
power-driven vessel, a fan blowing produces a similar low frequency content due to the 
blade rate of the fan as it rotates like a propeller. 
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2. Are Shallow Water Acoustic Signals Stationary? 

Like most random processes occurring in nature, the background ocean ambient 
noise is not stationary. In the example above, one noise source out of many may change 
its course or speed. So by one slight change in frequency, a change in the stochastic 
nature of the ambient noise will take place. However, as in many applications to linear 
filtering techniques, a process can be assumed stationary over a short period of time [5]. 

The benchmark signal for this thesis assumes that the signal is stationary for the 
entire 6 seconds. As the results will show (Figure II-2), assuming the data set is 
stationary over this time frame is sufficient to distinguish the transient signal from the 
background noise. Like the Benchmark Signal, the short-time stationarity assumption for 
underwater acoustics will apply. 



B. STOCHASTIC SIGNAL MODEL 

The first of two current methods to be discussed is an analysis in the time domain. 
We wish to transform the input signal through a filter to extract a transient from noise. 
To do this, we rely on the innovation process where the output of the filtered sequence is 
white noise. 

It is well known that a process whose complex spectral density function satisfies 
the Paley-Wiener Condition can be modeled as an Autoregressive-Moving Average 
(ARMA) process [5]. One of the properties of the signal model is that it can be reversed 
such that it becomes a whitening filter. Furthermore, since its transfer function H ca (z ) is 

minimum-phase, its inverse H ca ~\z) will be causal as well. Figure II-3 shows a diagram 
of the innovation representation. 
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Figure II-3 Innovation Representation of a Random Process, (a) Signal Model 

(b) Inverse Filter [5] 



1. Autoregressive - Moving Average (ARMA) Model 

Any regular stationary random process can be represented as the output of a linear 
shift-invariant filter driven by white noise [5]. Starting with the assumption that our 
input signal jr[n] in Figure II-3(b) is colored noise, we want to design a filter H(z) such 
that the result will be white noise. The innovation process will whiten the time domain 
signal except where the transients are located because it is not part of the innovation 
process. A second look at Figures II- 1 and II-2 will show this process. 

We use an autoregressive moving average (ARMA) model in order to represent a 
random process. The ARMA model involves both a nontrivial numerator and 
denominator and when compared to the AR or MA models alone, they often require 
considerably fewer parameters to match a desired power spectral density function [5]. 
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Without a priori knowledge of the autocorrelation function, we estimate it 
directly from the data. The Autoregressive model estimation can then be written as a 
weighted sum of prior values of the observations. Written as an equation, we get 



x[n] = -a x x[n-\]-a 2 x[n-2]- ,..-a p x[n- P ] , (2.1) 



where the error is 

p 

£[ri] = *[/?] - .?[«] = ^a k x[n — £] . (2.2) 

/t=o 

In equation 2.1, jc[n] is the estimate of the n' h element. As in many applications, 
we want to minimize the error by applying the Orthogonality Principle which states that a 
vector of weighting coefficients a can be chosen to minimize the mean-square error 

H{|,r — Jcj" } if the error is orthogonal to the observations. The final result is the least 
squares form of the Yule-Walker equations. 
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and S is the sum of the squared errors. It follows then that an estimate for the 
prediction error variance is 



a 







(2.4) 



where N s is the size of the data set. 

The observation matrix X is set up such that it assumes the observations 
outside of N s data points are zero, which is analogous to performing a rectangular 

window operation on the data set. Although this might not be the best method to use in 
terms of pole placement of the vector a, the matrix structure has a nice property since 
(X‘ r X) is Toeplitz. This property means that the matrix will be positive definite, which 
guarantees that the whitening filter is minimum phase and therefore is a causal stable 
system with a causal stable inverse. This property is important because we are interested 
in the innovation process. Finally, using the autocorrelation method of determining the 
Toeplitz correlation matrix (X* r X) , fast algorithms for determining the elements of a 
can be used such as the Levinson Recursion [5]. From the AR coefficients found from 
equation 2.4, the Moving Average coefficients b can then be computed by any number of 
algorithms such as the Basic Prony Method [5]. The ARMA routine used in this thesis 
uses an algorithm to determine an optimum number of poles and zeros, P and Q, 
respectively. 

As an example, take the Benchmark Signal and find the AR filter coefficients a 
and the MA coefficients b where 





■ i ' 




' 1 ' 


a - 




and b = 


b i 




. a p. 




_v 
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The causal ARMA filter is then 



^(z"‘) = 



b n +b,z 1 +b-,z 2 ■■■ + b n z Q 



-k u o n "l 



1 + ajZ" 1 + a,z" 2 H a p z P 



(2.5) 



To realize the time domain solution to the problem described in Figure II-3(a), we 



get 



x[n] = Yjb q \v[n - q] - £ - p] . 

< 7=0 p = 1 



( 2 . 6 ) 



Inverting the model as in Figure II-3(b) the causal inverse ARMA filter is simply 



h arm A^ = 



1 + a,z 1 +a 2 z 2 H h a p z P 

b 0 +b y z~ l +b 2 z~ 2 +--- + b Q z~ Q 



(2.7) 



The time domain realization to the innovation process is finally. 



w[n] = Yjd p x[n - p] -^\b q w[n - q ] . 

P=1 <7=0 



( 2 . 8 ) 



Running the Benchmark Signal through the inverse ARMA filter (as depicted in 
the model in Figure II-l) results in the innovation process, which is white except where 
the transients occur in time. This is because the ARMA model seeks only to model the 
colored noise as the innovation process. Since the colored noise input has a stationary 
localized frequency content, the ARMA model can do a good job at placing poles to 
model the signal. All other signal content that is statistically dissimilar to the colored 
noise is rejected and appears as a “spike”. 

As a check, in listening to the fan blow, a low frequency “hum” is distinguishable 
and the 1 1 tapping sounds are also recognizable. After the innovation process, the 
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background noise sounds like a heavy rain hitting the ground. When plotted the 
transients, indistinguishable in the time domain before the innovation process, are now 
quite visible against the white noise. Figure II-4 shows the Benchmark Signal before and 
after the innovation process. 
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2. Strengths of the ARMA Model 

Since the innovation process is based on a statistical model, the frequency content 
of the signal does not affect its detection, even if the signal contains the same frequencies 
as the noise. The ARMA model provides a representation of the noise and it can detect 
signals that are not compatible with the model. The ARMA method picked up on all the 
transients and the height of the spikes of the output signal is proportional to the 
“loudness” of the audible signal, which might be an attractive feature to compare the 
magnitude of two signals. This technique seems to perform well on very short transients 
and it can resolve signals close to each other. Notice, for example, the cluster of signals 
just before the third second in Figure II-4 (b), where the ARMA model was actually able 
to distinguish between three transients less than one-tenth of a second apart from each 
other. 



3. Weaknesses of the ARMA Model 

The reverse ARMA (innovation) process is a procedure that produces white noise 
from a colored noise input. This process, as we shall see, is very effective in detecting 
transient signals that are not part of the white noise model, appearing as spikes above the 
noise threshold as seen in Figure II-4 (b). However this method, although performing 
well on both narrowband and broadband signals, has its limitations when the power of the 
transient signal is at or below the modeled noise power. As we shall see, the fundamental 
idea of the methods presented in this thesis is to lower the white noise power of the 
innovation process in order to extract transients with smaller signal power. 
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c. 



TIME-FREQUENCY ANALYSIS 



The second method to be discussed introduces a technique that combines both 
time and frequency analysis. The process is like a sliding window in time taking a 
snapshot of the frequency content of a signal, assuming that the statistical properties do 
not change within the time duration. To view the results, frequency is plotted as a 
function of time and is called the spectrogram. For USW, a similar technique is used 
called the “lofargram,” which is currently the standard method for sonar signal 
processing [5]. 



1. Spectrogram 

The Spectrogram is based on the discrete-time Short-Time Fourier Transform, 
which is an extension to the basic Discrete-Time Fourier Transform and is defined as 



X(n,co) = ^ x(m)w(n-m)e Jwm , 

m =~<> ° 



(2.9) 



where 



w (/2 - m) is the sliding window function and 



X(n,0))= £ x(m)e~ iam is the DTFT. 



Completely analogous to the DTFT, the discrete-time STFT is continuous in 
frequency so for digital processing, a method to discretize frequency content like the FFT 
can be implemented to make this method computationally practical. Figure II-5 shows a 
diagram of the STFT process. Along the horizontal axis is a sliding window process as a 
function of time. At each T lasl + A? , a “snapshot” of the frequency content of the signal 

at that instant in time is taken. Another way to visualize the process is to look at the 
frequency domain. Along the vertical axis is a set of filter banks, each of equal 
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bandwidth A/ . The signal is passed through the filter bank at each r last + At and the 

spectral content within the particular bandpass filter is calculated. The latter pictorial 
representation of the STFT is particularly useful in drawing an analogy to the Continuous 
Wavelet Transform discussed next. 
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Figure II-5 Time-Frequency Plane Corresponding to the Short-Time Fourier 

Transform [7] 
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The (discrete) STFT on the Benchmark signal is provided in Figure II-6. The 
darker colors on the spectrogram represent frequencies with high power. Notice the low 
frequency content below 300 Hz is dominant throughout the entire 6-second sample. 
Constant frequencies at 600 Hz and 900 Hz are also noticeable. The most revealing 
information on the spectrogram is the vertical lines in the latter half of the signal 
indicating the presence of transients. What’s more, the discrete STFT can provide 
information on the type of transient it is. In this case, the transients are broadband 
signals. In addition to the standard Benchmark Signal, three narrowband transients were 
added to the first second to study the effects of narrowband transient resolution. The 
spectrogram located two of the three signals, one at around 5 kHz ± 500 Hz and one 
around 2 kHz ± 500 kHz. The third narrowband transient, however, was not detected 
because it was masked by the low frequency noise. 
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Figure II-6 Spectrogram of Benchmark Signal 



2. Strengths of the Spectrogram 

Among the many strengths of the spectrogram is its colorful visual representation. 
On a single plot, it provides information in both the time and frequency domains. This 
representation therefore gives not only the time at which the transient signal occurred, but 
also its frequency content, which may be vital in classifying the source. For applications 
in sonar signal processing, the spectrogram is commonly referred to as the “lofargram”, 
which provides a useful “fingerprint” mapping of unique acoustic sources [5]. From this 
display alone, useful information such as detecting dominant narrowband tonals. 
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harmonics and blade rates can be seen. Seeing the frequency content of a contact change 
over time can also indicate a dynamic event such as a change in course or speed. 

The spectrogram works best on broadband transient signals because there is a 
larger contrast of darker colors in a vertical line than there is in the single point 
representation of a narrowband signal. In addition, we have shown that it is possible to 
have a narrowband transient signal masked by the frequency content of the background 
noise. Unless the background noise is completely white, at least some of the frequency 
content of a broadband noise is sure to appear. 

Compared to the ARMA model, the STFT only assumes stationaritv over the 
duration of the sliding window only. Recall that the ARMA process attempts to model 
the signal through the statistics of the entire 6-second data set. If it were to assume a 
shorter data set to compute the autocorrelation matrix, the whitening model will be 
different and may or may not provide better results. A further discussion on the window 
size will be discussed in the final chapter. 

3. Weaknesses of the Spectrogram 

Like the innovation process, the discrete STFT must assume stationarity for the 
duration of the windowing function. There exist natural and man-made signals for which 
its spectral content is changing so rapidly that finding a proper windowing function is 
difficult because there may not be any obtainable time interval for which the signal is 
stationary [6]. The result is a transient signal may go completely unnoticed. 

Even if it is a good assumption that a short time duration is stationary, there is a 
tradeoff to the resolution in the frequency domain and vice versa. This is called the 
Uncertainty Principle or Heisenberg’s Inequality, which states that resolution in time and 
frequency is bound by the inequality 



An 



(2.10) 



with At and A / the duration of the signal in time and frequency respectively. 
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If localization in time is necessary, At can be made small but according to the 
inequality (2. 10) A/ will become greater and therefore the spectrogram will lose 
frequency resolution proportionately. 

To further illustrate this principle, the two visible narrowband transient signals in 
Figure II-5 are about one-tenth of a second in duration and contain only one frequency 
each, 5 kHz and 2 kHz respectively. From the figure, both transients are stretched 
vertically, blurring the actual frequency spectrum while the time duration at which they 
occur remains short. If a broad window were used in the discrete STFT, the frequency 
would resolve to a single point while the time of the signal would be stretched and 
distorted. 

The spectrogram merely shows the frequency of a signal during a time duration. 
Since the spectrogram is not based on statistical modeling, it cannot “extract” a signal 
from the noise and therefore is unable to distinguish between signal and noise if they both 
have the same frequency content. This is particularly important in cases where an 
acoustic transient signal is masked by the colored noise. Consider the Benchmark Signal 
spectrogram in Figure II-5. Also included in the signal is a sine wave at a frequency of 
100 Hz. Since the power of the colored noise is highly concentrated at frequencies below 
300 Hz, the transient signal goes undetected. Unmasking transient signals in this region 
become important for sonar operations because much underwater activity may take place. 

Finally, while the innovation process forces a signal into a white noise process in 
order to detect transients, the spectrogram takes advantage of the signal’s frequency 
content over time. If the input signal is white noise, nothing is distinguishable since the 
frequency content is spread evenly over all frequencies for all time. 
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III. WAVELET ANALYSIS 



Since we know that the STFT is resolution-limited in both time and frequency 
simultaneously, we wish to investigate other means to localize transient signals in time 
while somehow still being bounded by the uncertainty principle. In this section, we start 
with a discussion of the Continuous Wavelet Transform (CWT) as it relates to the Short- 
Time Fourier Transform. From there, we discuss the development of the Discrete 
Wavelet Transform, its usefulness with filter banks, and the idea of multiresolution 
analysis, which is the basis of the methods of detecting transient signals. 

A. CONTINUOUS WAVELET TRANSFORM 

To describe the Continuous Wavelet Transform, w'e start with the equation for the 
Short-Time Fourier Transform 



The integral can be viewed as an inner product operating on data inside the 
windowing function, which measures the “similarity” between the function x(t ) and a 
sinusoidal basis function. The result turns out to be the energy level at each specific 
frequency. STFT analysis is useful for narrowband periodic signals because of its ability 
to detect strong similarity between the input signal and the sinusoidal basis function at a 
particular frequency, but it has its limitations to the larger family of non-stationary, 
broadband signals. Based on Fourier Transform restrictions, it is desirable to develop a 
different transform that can overcome the resolution limitation. 

In order to investigate a different set of transforms that exhibit these properties, 
w'e need to look at other possible basis functions. To begin, assume that the width of the 




( 3 . 1 ) 
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frequency resolution of the transform A/ is no longer constant but rather the ratio of A f to 
the frequency / , resulting in 



/ 



= c. 



(3.2) 



This constant ratio means that as the analysis of frequency content of a signal 
increases, so does the width of the frequency resolution. Looking back at the diagram of 
the time-frequency plane of the STFT (Figure II-5), the vertical axis represents a set of 
constant-width bandpass filters in a filter bank. In the case of the CWT, the constant- 
width frequency resolution is replaced by a filter bank of constant relative bandwidth 
(called “constant-Q” analysis). This process, in essence, allows the resolution in At and 
A/ to vary in the time -frequency plane in order to obtain a multiresolution analysis, 
which will be discussed in greater detail in the next section. Another way to look at the 
filtering process is that the filter bank used is not linearly spaced along the frequency axis 
(as for the STFT case) but spaced evenly over a logarithmic scale [7]. Figure III. 1 shows 
this process. Now substituting (3.2) into the Heisenburg Inequality (2.10) gives 



At = 



1 

4 7tfc 



(3.3) 



By looking at equations (3.2) and (3.3), one can see that as the frequency becomes 
greater, the resolution in time decreases and the resolution in frequency increases. 
Conversely, as the frequency becomes smaller, the resolution in time becomes larger and 
the resolution in frequency becomes smaller. This is a fundamental property of the CWT. 
The time resolution of the CWT becomes arbitrarily good at high frequencies while the 
frequency resolution becomes arbitrarily good at low frequencies. 
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Figure III-l Filter Division in the Frequency Domain, (a) Uniform Coverage (STFT) 

(b) Logarithmic Coverage (CVVT) [7] 



The definition of the Continuous Wavelet Transform is thus 
CWT X (r, a) = J x(t)h\,r ( t)dt , 



where 



h\At) 




( t-r\ 



a j 



(3-4) 



(3.5) 



23 



is the basis function called a wavelet. 

A description of a wavelet function will be addressed in the next section. For 
now, some of the properties of the transform will be discussed. From the definition, it 

can be seen that a is the “scaling” factor described in equation (3.2) and 1/ ^J\a\ is used 

for energy normalization. The numerator of the wavelet function shows the time - shift 
property of the sliding window function. 

By way of analogy, the STFT is plotted as frequency vs. time as a spectrogram. 
The CWT is a function of time and scale and therefore will be plotted as scale vs. time as 
a “scalogram”. Figure III- 1 is a comparison of a CWT to a STFT of the same signal. 
Notice that the CWT of the delta function converges to a single point at t = t 0 and the 

energy is spread out with increasing scale. The STFT of the delta function shows the 
resolution limitation in the time domain with equal energy at all frequencies. For the 
sinusoidal signal of three frequencies at f 0 , 2 / 0 and 4 f 0 , the STFT shows the constant 

frequency in time, bound by the resolution limitation in the frequency domain. For the 
CWT, however, we see the effect of the logarithmic filter bank where the low frequency 
resolution is narrowest and continues to double in size as the frequency increases. Figure 
III-2 shows clearly the principle that the time resolution of the CWT becomes arbitrarily 
good at high scale while the frequency resolution becomes arbitrarily good at low scale. 
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Figure III-2 Comparison of the CWT to the STFT [7] 
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B. DISCRETE WAVELET TRANSFORM 



At this point a discussion of the properties of the Continuous Wavelet Transform 
has been made as it relates to the STFT and the benefits in localization in both time and 
scale have been presented. The CWT also reduces the number of coefficients necessary 
to define the signal. During the discussion, ideas of orthogonality, energy distribution, 
logarithmic filter banks, scaling functions, and wavelet basis functions emerged as 
important properties of the CWT. This section will continue to build on these concepts as 
they apply to the Discrete Wavelet Transform (DWT). 

For reasons of clear analysis and better understanding, we begin with the 
definition of a linear decomposition. A real- valued function can be expressed as a linear 
decomposition such that 



/(0 = 2<mM0- (3.6) 

l 

Parameter a i is the real-valued expansion coefficients, and y/(t) is a set of real- 
valued functions of t called the expansion set. When the expansion is unique, the set of 
functions y/ t (t) is called a basis and it is used to represent fit). If the basis is 
orthogonal, i.e. 



(y/ k ,y/,) = ji// k (t)il/,(t)dt = 0, k*l (3.7) 



then the coefficients can be determined by the inner product 

a k =(fit),vr k ) = \fit)i/r k (t)dt. (3.8) 



For the wavelet expansion, however, a two-parameter system is used such that the 
basis used to represent fit) is now 
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(3.9) 



f(0 = 'Z'Za Jk y/ jk (t), 

k j 



where 



j represents the integer scaling index of the scaling function, 
k represents the integer time translation, 

if/ j k (?) are the wavelet expansion functions that usually form an, 
orthogonal basis, and 

a j k are the set of expansion coefficients. 



The wavelet basis set y/ jk (t) can be defined in terms of scaled and translated 
versions of the Mother Wavelet y/(t) , i.e. 



¥ j,k (0 = 2^- y/(2 J t - k) . 



(3.10) 



Here, the term 2 2 is a normalizing factor as j increases. This property described 
in equation (3.10) is called the multiresolution condition, which states that if a set of 
signals can be represented by a weighted sum of <p(t - k) , then a larger set (including the 
original) can be represented by a weighted sum of (p{2t-k). In other words, “if the 
basic expansion signals are made half as wide and translated in steps half as wide, they 
will represent a larger class of signals exactly or give a better approximation of any 
signal” [9]. 
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1. Signal Spaces and Multiresolution Analysis 

To discuss the DWT, it is best to start with the idea of resolution to define the 
effects of changing scale. To that end, we start with a scaling function (pit ), which then 
will define the wavelet basis function \f/{t ) . We define a set of scaling functions in terms 
of integer translates of the basic scaling function by 

<p k (t) = (p{t-k) keZ and (peL 2 . (3.11) 

where the L 2 space is defined as the set of all functions / ( t ) with finite energy, 

i.e., 



j\f(t)\~dt <°°. 

— OO 

The subspace spanned by the functions (p k ( t ) is v 0 such that 

v’ 0 = Span{(p k (0, keZ}. (3.12) 

That is to say, the subspace v 0 is spanned by the closed set of all functions <p k ( t ) 
over all integer values k . From equation (3.12) we can say that a function fit ) within 
the subspace v 0 can be expressed as a linear combination of the set of functions (p k (r) 
such that 

f(t) = Yj a k<Pk( t ) for any f(t)ev 0 . (3.13) 

k 

From here, we consider a two-dimensional family of functions and derive a similar 
linear combination as in equation (3.13). We start with the definition of a scaling 
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function parameterized by scaled and translated versions of cp(t) similar to equation 
(3.10), resulting in 



<p jk (t) = 2'-<p(2 J t-k), keZ (3.14) 

where 

v ] = Span{<p jk (/), k € z\. (3.15) 

In the same way v 0 is expressed as a linear combination of the set of functions 
(p k {t) in equation (3.13), f(t)e can be expressed as 

f(t) = 'Za k <p(2 j t + k). (3.16) 

k 

From here we turn to a useful property of the wavelet expansion set [9], There is 
a basic constraint of multiresolution analysis, which requires a “nesting" of the spanned 
spaces as 



. . . c v_ 2 c v_j cv 0 c V', c v, c . . . L 2 , 



or simply 



Vj cv ;+1 for all jeZ. (3.17) 

By this observation, in order to ensure that the elements in Vj are scaled versions 
of the next higher space v j+] , we add the constraint. 
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/(Oe Vj 



<=> f(2t)ev 



(3.18) 



j+i • 

This is done to ensure that analysis and synthesis can be performed. That is to 
say, if we conduct a wavelet decomposition at a certain resolution (or scale), there exists 
a reversing process whereby a signal can be reconstructed (synthesized) to its original 
state. Figure III-3 shows a Venn diagram of the nested vector space property of 
multiresolution analysis. 




This set of nested vector spaces now meets the multiresolution condition, which 
states that if a set of signals can be represented by a weighted sum of (p ]k {t), then a 

larger set can be represented by a weighted sum of (p j+Ue (t) ■ In order to relax notation, 
we drop k and assume that (p,(t) is dependent on the time translation. The nesting of the 
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spans as shown in Figure III-3 and by equations (3.17) and (3.18) is obtained by requiring 
that <p 0 (t ) is not only in the space v 0 but also in v,, which is spanned by the function 

<Pi (?) . From equation (3.16), we can see that that the function <p 0 (t ) within the subspace 
v 0 can be expressed as a linear combination of the set of functions <p { (?) such that 

^o<» = X a n^i(0- (3.19) 



Substituting <p 0 (t) and (p x {t) using equation (3.14) we obtain an equation in terms 
of scaled and translated versions of the scaling function as 

<p(t) = ^ h(n)y[2<p(2t — n), ne Z . (3.20) 

n 



This equation shows a telescoping structure of the scaling function and is called 
the Multiresolution Analysis (MRA) equation, where h(n) represents the weighted 
coefficients called the scaling function coefficients. 

We now have a recursive scaling function spanning a subspace, which is also a 
subset of the next higher scale. There is one more feature of a signal that does not use a 
scaling function, but rather uses a different set of functions, called y/(t). This set of 
functions is defined as a basis set that spans the differences between the spaces spanned 
by the various scales of the scaling function [9]. As we shall see, these are the set of 
wavelet functions. 

Although wavelets and their corresponding scaling functions are not required to 
be orthogonal, for the purposes of this thesis, we will limit our discussion to the set of 
orthogonal basis functions. The advantages of having an orthogonal basis are that the 
wavelet expansion coefficients are easy to calculate and are governed by Parseval's 
theorem which allows the partitioning of signal energy in the wavelet transform domain, 
like that of the classic Fourier transform domain. In vector function space notation, we 
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define w as the orthogonal complement of v ; in v j+1 . That is to say, } is the set of 
all functions that span v j+l and are not included in v, . v ■ is also be considered 
uncorrelated with w ; and the union of the two vector spaces span v +1 . Mathematically, 
all functions spanned by v j are orthogonal to vw . Therefore the basis functions (p ] k and 
y/ ,j are orthogonal as 



(<Pj.t (0, Wj, ( 0 ) = J (pj. k (0 Vjjt (0 dt = 0 . (3.21) 

Starting at an arbitrary scale spanned by the scaling functions in v. , say y =0, 
we can write 



v 0 c v, c v 2 c • • • c Is . 

Now from the orthogonal relationship between v 0 and w 0 , we get 

v,=v 0 ©w 0 , (3.22) 

where © indicates a union operator. 

Finally, due to the “telescoping” nature of the subspaces v j , j e Z and the 
orthogonality of the corresponding Wj , je Z , all of L 2 can be expressed as 

L~ = v 0 © u' 0 © vv’j © u’ 2 © • • • . (3.23) 

Figure III-4 shows the Venn diagram of the entire vector space spanned by L 2 . 
Assuming that our arbitrary starting point in the infinitely telescoping set of nested vector 
spaces spanned by the scaling functions is r , we can see that w . is orthogonal to v } 
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and a subset of v j+l as shown in Figure III-5. These wavelets reside in the space spanned 

by the next narrower scaling function and therefore the multiresolution principle applies. 
Since w 0 c v,, the set of wavelets in w 0 can be represented as a weighted sum of the 
next higher scale <p(2t ) resulting in 



where /i, (;j) represents the weighted coefficients that define the wavelet. 

It can be seen that in equation (3.24) the wavelet function is expressed in terms of 
a linear combination of scaled and translated versions of the scaling function, which 
implies a relationship between the scaling function and the wavelet function. Indeed, by 
taking equations (3.20) and (3.24), it can be shown that h(n) and /i, ( n ) are related by [9] 



neZ, 



(3.24) 



n 



hi(n) = (-1)" h(-n ) . 



(3.25) 
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Figure III-4 Scaling Function and Wavelet Vector Spaces [8] 
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2. Implementation of the Discrete Wavelet Transform 

We now have defined a wavelet system in terms of its scaling function and 
Mother Wavelet. As it turns out, it is not necessary to actually know <p(r)and y/(t ) . All 
we need are the filter coefficients to the MRA equation in order to perform a discrete 
wavelet decomposition. This section introduces the Discrete Wavelet Transform (DWT) 
and how it can be implemented. 

In the last section we described the whole vector function space as a union of the 
coarsest scale with all wavelet subspaces such that 

L~ = v 0 © w’ 0 © Wj © w 2 © • • • . 



Based on this fact, we can now express a function g(t) as a series expansion in 
terms of the scaling function <p k (t) and the wavelet functions y/ ]k (t) giving the general 
form of the equation 



g(0 = E c ;o(^)^yo,i (0 + X Yj d • ( 3 - 26 ) 

k j=j0 

Here, jO is the arbitrary position for the coarsest scale whose space is spanned by 
the scaling function (p Q k (t ) . Notice that the inner summation on the second term begins 
at jO , indicating that the wavelet coefficients begin at the same level as the current scale 
level. If the conditions of orthogonality are met, these wavelet coefficients can 
completely describe the function g(t) and are useful for analysis, denoising, compression 
and perfect reconstruction of the original signal [10]. Assuming that the wavelet system 
is orthogonal, the coefficients c .(&) and dj(k) can be found by the inner product of (p j k 

and iff jk respectively. 

Performing the inner product to find c } {k) we have 
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(3.27) 



Cj (k) = (g(0, (p>' k (0) = J g(t)<Pj, k (0 dt 
= {g(0 2*<p(2 j t-k)dt. 

Restating the MRA equation, 

(Pit) = Yj h(ny/2q>(2t - n ) , 

n 

we substitute in the next higher scale and translate the function by k, resulting in 
(p{2 1 t-k) = YHn)j2(p(2(2 j t - 2k - n)) . 

n 

By making a variable change m-2k + n, we get 

(p(2’t - k ) = Y h ( m ~ 2k)^2<p(2 J+l t-k ) . (3.28) 

m 

Now substituting (3.28) into (3.27) it can be shown that 

c y (*) = 5>(m-2A')fg(02 2 <p(2 i+l t - k)dt . (3.29) 

m 

Finally, the integral expression is simply the inner product 

(s(0,^ + i. m (0) = c y+1 (m ) , 

and therefore 
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(3.30) 



c j (*) = 2 - -k) c j+i ("0 • 

A similar derivation for finding the coefficients d ( k ) yields 

dj (k) = /z, (m - 2 k)c j+l (m) . (3.31) 

m 

Equations (3.30) and (3.31) show that the scaling and wavelet coefficients at the 
next lower scale j can be computed by convolving the scaling coefficients at the current 

scale j + l with the time-reversed weighting coefficients h(-n) and fyf-n), then 
downsampling the result. Figure III-5 shows a diagram of this process. 




Figure III-6 Two-Stage, Two-Band Analysis Tree [8] 
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3. Analysis in the Frequency Domain 

In the following section, it will be shown that the coefficients of the Haar Wavelet 
filter coefficients h 0 (/?) are a lowpass filtering operation. Likewise the coefficients fyO?) 

are a highpass operation [9]. Together, they form a Quadrature Mirror Filter Bank 
(QMF) system. In addition, a reverse filtering operation can be performed allowing for 
perfect reconstruction of the original signal, also known as a synthesis procedure. This 
makes sense because for any transform, a reverse operation must be obtainable for the 
transform to be widely used. 

In the frequency domain, the signal coming in is processed through a filter bank 
that splits the frequency content into two parts. The high frequency band of the signal is 
represented by the wavelet coefficients also known as the details. The low frequency 
portion of the signal is represented by the scaling coefficients also known as the 
approximations. The lowpass portion can then be processed through the same filter bank 
again. Theoretically, lowpass filtering the scaling coefficients can be repeated 
indefinitely, but in practice is limited to the number of data points in the signal. This 
limitation exists because every time the signal is passed through the filter bank it gets 
decimated by a factor of 2, which effectively means the number of data points in the 
signal is cut in half. 

The first stage splits the frequency band in half, the second stage splits the low 
frequency portion into quarters, and so on. The resulting “analysis tree” makes up a 
logarithmic filter bank as seen in the CWT. A diagram of this process can be seen in 
Figure III-7 as it compares to a constant width filter bank. 
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Figure III-7 (a) Constant Width Filter Bank (b) Frequency Bands for the Analysis 

Tree [8] 



39 




C. HAAR WAVELET 

Up to this point we have studied the Wavelet Transform as a linear 
decomposition. We now move from the abstract to the application. Equations 3.30, 3.31 
and Figure I1I-6 all describe the process of decomposing a signal into successive details 
and approximations. Now looking at the decomposition as an application to signal 
processing, we see that it is a series of simple filtering and downsampling operations. 
For our particular application, the actual numerical values of the decomposition 
coefficients are not important, only their relative values and so the coefficients can be 
normalized. By viewing the DWT as a simple signal processing procedure, we have 
essentially found a way to split and downsample a signal into mutually independent 
signals. As we shall see, this procedure provides useful properties in detecting transient 
signals by removing background noise. 

1. Haar Wavelet Parameters 

We are now ready to discuss the actual wavelet filter coefficients h 0 (n ) and 
. There are many wavelet basis sets to choose from and a discussion on choosing 
specific types will be discussed in the final chapter. For the purposes of this thesis, we 
will choose the Haar Wavelet system and derive all transient detection methods and 
conclusions based on this set. Admittedly, this is the simplest of all wavelet basis sets but 
the results will show a marked improvement over current detection schemes. In addition, 
it is more intuitive to describe the filtering processes in the upcoming chapters using Haar 
wavelet coefficients. Once the transient detection methods have been discussed, a more 
general case can be easily substituted. 

The Haar Wavelet coefficients are defined as 




(3.32) 
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for the lowpass filter and 



* I (n) = {A 1 (0)A(l)} = 



1 

V5 



-l 




(3.33) 



for the highpass filter. 

A more detailed discussion explaining how the Haar Wavelet meets the necessary 
conditions of a Wavelet system is presented in Appendix A. To find the shape of the 
scaling function, recall the MRA equation (3.20) 

(pit) = ^h(n)y[2<p(2t ~ n) . 



Substituting the scaling function coefficients, the scaling function becomes 

(pit) = (pilt) + (pill -n). (3.34) 

Equation (3.34) clearly shows that the Haar Scaling function is a unit-width, unit- 
height pulse function and is equal to the sum of two scaled and shifted versions of itself. 
To find the shape of the Haar Wavelet, we use the wavelet equation (3.24) 

V(t) = 'ZK(n)S<p(2t-n). 



Substituting the Haar Wavelet function coefficients, the above equation becomes 

y/it) = (pi2t)-(pi2t-\). (3.35) 



In the same manner as the scaling function, the shape of the Haar Wavelet is the 
sum of two scaled and translated versions of (pit ) , except now the coefficient on the 
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second term is negative reflecting the negative coefficient in (n) . A diagram of the 
Haar Scaling and Wavelet functions is shown in Figure III-8. 




Figure III-8 (a) Haar Scaling Function (b) Haar Wavelet Function 
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For the analysis of both transient detection procedures, it is vital to understand the 
effect of the wavelet decomposition filtering process in the frequency domain. We start 
out by viewing the Haar Scaling filter coefficients in the z - domain as 

K(z~') = ~r( 1 + z_1 )- (3.36) 

V2 

Substituting e J " for z we get the frequency response of the FIR filter 




(3.37) 



Likewise, the frequency response of the Haar Wavelet filter coefficients is 



H x {co) 





(3.38) 



It is important to note here that the frequency response of the scaling coefficients 
h(n ) (also known as h 0 (n )) is that of a lowpass filter with its cutoff frequency at the 
Nyquist frequency rate ( H{n ) = 0). It can also be shown that the frequency response of 
the wavelet coefficients h x {n) is that of a highpass filter with its maximum at 

H(n) = V 2 and its zero magnitude at H( 0) . This is the property that allows the input 
signal to be filtered into a “constant-Q” filter bank resulting in equal highpass and 
lowpass portions preceding every downsample as described in Figure III-7. Figure III-9 
shows the frequency response of the PR filter bank used for the Haar Wavelet system. 
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Figure III-9 Frequency response of Haar Scaling and Wavelet Functions 



2. Haar Wavelet Summary 

For the DWT, the Haar Wavelet and Scaling functions themselves are only 
important in terms of describing the wavelet decomposition process heuristically. The 
only parameters needed to apply wavelet analysis are the scaling coefficients to the 
Multiresolution Analysis equation. Table III-l shows a summary of the Haar Wavelet 
system. 

Indeed, the scaling coefficients h 0 (n) and the wavelet coefficients /*,(») work 
together to form a perfect reconstruction filter bank. When a signal is processed through 
the FIR filter /^(n), it is being passed through a highpass filter. When the signal is 
downsampled, the set of values is the wavelet decomposition coefficients called the 
“details”. Likewise when the same signal passes through the FIR filter h 0 (n ) , it is being 

processed through a lowpass filter. After the signal is downsampled, the set of values is 
the scaling decomposition coefficients called “approximations”. The approximations can 



44 



then be decomposed again using the same filter bank and downsampling indefinitely 
(theoretically) to allow for infinite (theoretically) resolution in both time and scale 
(Figure III-6). 

Finally, all details at different scale and the approximation at the smallest scale 
are orthogonal to each other. This property, as we shall see, implies that all signals at the 
output of all filter bank-downsample phases are uncorrelated to each other and is the 
basis for the first transient detection method. 




Table II-l Haar Wavelet System Summary [15] 
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IV. METHOD 1: INNOVATION - DWT PROCESS 



With the desirable properties of the DWT, we focus on the first method to 
enhance transient detection. This section will present a signal detection design based on 
decomposing the innovation process through the Haar Wavelet system. The desired 
result will be to remove noise through the filtering process while keeping the signal 
intact. 



A. METHOD 1 DERIVATION 

As mentioned in the previous section, the first method takes advantage of the 
orthogonality of subspaces in wavelet decomposition. To say that the details and 
approximations at each stage are orthogonal implies that they are uncorrelated to one 
another. As a proof, consider the whitened input signal processed through one stage of 
the Haar Wavelet system as in Figure IV-1. Note that e LF [/] and e HP [l] are using the 
same index l because they are produced from the same input and are decimated by the 
same value. 
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Figure IV-1 Wavelet Decomposition on Whitened Signal 



From the previous figure, we can express the approximations e LP [l ] and the 
details e HP [l] in terms of the input signal e[n] such that 



e LP [l] = ^{e[2l] + e[2l-l]) and 


(4.1) 


e HP [l] = \{e[2l]-e[2l-l]). 


(4.2) 



To determine the correlation between the two signals, we perform the cross 
correlation operation, i.e. 

E{e LP [ l ] ■ e HP [/]} = i E{(e[2l] + e[2l - 1]) • {e[2 1] - e[2 1 - 1])} (4.3) 
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Carrying through the multiplication, the expression yields 



E{e LP [l] -e HP [l ]} : = ±e{ \e[2lf }- ± e{ \e[2l - 1]| : 2 }= 0 



(4.4) 



Since the signal e[n] is wide sense stationary, £ , jje[2/]| j is equal to 

£{|e[2/ - 1]|‘ } and equation (4.4) is zero and hence, e LP [l ] and e HP [l] are uncorrelated. 

The idea of Method 1 is to perform a DWT on the whitened signal. This setup 
passes a white noise process through the wavelet decomposition where we expect white 
noise at each resolution stage. Based on the properties of the Discrete Wavelet 
Transform and the proof at the beginning of the chapter, all the details and the 
approximations at the same multiresolution stage are uncorrelated. The description here 
is a denoising process where we attempt to lower the noise power by successive filtering 
and downsampling procedures. 



B. METHOD 1 SUMMARY 

Method 1 uses the Haar Scaling and Wavelet filter coefficients for the 
multiresolution filter bank system and is shown in Figure IV-2. Keep in mind this 
method assumes that the DWT is operating on a white noise signal, which means that the 
use of the ARMA model is used prior to the first stage of the filter bank. If the 
innovation process detects a transient signal, note that there is no need for this algorithm. 
What this method attempts to do is go beyond the detection capabilities of the innovation 
process by removing excess noise in order to enhance the robustness of the detection 
process. 
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V. 



METHOD 2: DYVT - INNOVATION PROCESS 



We now take a different approach to the same problem using a derivative of the 
Haar Wavelet and its multiresolution capability. As shown in Figure IV-1, Method 1 
took the wavelet decomposition after the innovation process. Now we will investigate a 
method where we perform the multiresolution filter bank process before the innovation 
process as in Figure V-l. The idea of this method is to investigate the effectiveness of 
another denoising scheme. 

Method 2 takes a signal and processes it through the wavelet decomposition. By 
the properties of the DWT, the signals at each channel will be mutually uncorrelated. At 
this point, the decomposed signal is still not whitened. In order to whiten the input, there 
must be an independent ARMA filter for each channel because the same set of ARMA 
filter coefficients is only useful for one channel. Fortunately, it turns out that just as there 
is a relationship between the filters in the wavelet filter bank, there is a corresponding 
relationship between ARMA filters in the innovation filter bank and is dependent on the 
Haar Wavelet filter coefficients. 

This chapter will discuss a procedure to determine a set of whitening filter 
coefficients from the original set of ARMA filter coefficients. The first section starts out 
with modified Haar Wavelet filter bank and begins the procedure to determine all the 
ARMA filter coefficients. 
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A. 



METHOD 2 DERIVATION 



Figure V-l shows a block diagram of the multiresolution innovation process. 
First, note that the multiresolution filter bank is a generic highpass/lowpass filter pair that 
are orthogonal to one another. Instead of the Haar Wavelet, the filter coefficients are 

multiplied by a scalar value of — so that 



;i °(”) = 7 o k °' 



Haar 



( n) = 



9 9 



(5.1) 



and 



Mn) 



^ th.Haar ( n ) — 



9 



-1 

2 



(5.2) 



The new set of highpass and lowpass filters can be viewed as an averaging 
(lowpass) filter and a difference averaging (highpass) filter. The innovation at the output 
of the whitening filter bank are orthogonal to each other where the “tildes” represent the 
highpass filter channels, or in terms of the DWT, the details. The subscripts represent the 
relative sampling frequency in the sense that e k is decimated by a factor of 2 k . The 

dashed lines in Figure V-l represent a set of cascaded lowpass filters that are 
intermediate steps for further decomposition, yet can also provide useful information. 
Although the decomposition can go on indefinitely, there comes a point where a short 
transient signal is averaged excessively with the background noise. By experimentation, 
three levels of resolution were found to be the most useful in this thesis. 
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It is also important to note here that the wavelet and ARM A filter coefficients are 
interrelated, and the model is designed exclusively for the Haar wavelet and its 
derivative. A more general form of Method 2 is possible, allowing the use of any wavelet 
system, but is a subject for further research. 

Consider the ARMA process without multiresolution as in chapter 2. For notation 
simplification, we define the ARMA filter as a transfer function operator on a white noise 
process such that. 



y 0 [n]^^rr:e[n], (5.3) 

Mz ) 

where 



A(z ‘) = 1 + a { z '+ — 1 -a p z P , 


(5.4) 


B(z-') = b 0 +b lZ - 1 +--- + b Q z- Q , 


(5.5) 



and e[n ] is white noise with the second moment being 



E{\e[nf}=(r w 2 . (5.6) 

From this definition we define y,[/] as the output of the lowpass/downsample 
procedure, also known as the approximation of the DWT. The goal of this section is to 
determine the ARMA filter coefficients of y,[/] as they relate to the ARMA coefficients 
of the original signal. Figure V-2 shows a simplified diagram of Method 2 where we 
consider a single stage of the approximation channel. Once the procedure to determine 
the ARMA filter coefficients at this channel is complete, we will go back and consider 
the ARMA model for the detail channel. 
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Looking again at Figure V-2, we want to expand the innovation process in the 
sense that we want to decompose a signal y 0 [«] as two independent signals as in Figure 

V-2. The transfer function f[z~ ] ) is the lowpass averaging filter described in equation 
(5.1) so that 



+ (5.7) 

The sequences e 0 [n] and e,[/] in figure V-2 are white noise processes where 
£,[/] is at half the sampling rate as e 0 [n] . Note here that we define different white noise 
processes throughout the derivation and the terminology can be confusing. The white 
noise signal e,[/] is the final state of the multiresolution filtering process. At this point 
there is no guarantee that e 0 [n ] and e x [l ] are uncorrelated and is subject for future 

research. However, it will be shown that the method does provide expected 
improvements in denoising. 
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Figure V-2 Single Stage of Multiresolution Innovation Process 
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Starting with the ARMA model of the original signal, we split the original signal 
into even and odd parts. We will define the equation now followed by the proof. The 
goal for this setup is to produce two white noise processes that are independent of one 
another. From the ARMA model in equation (5.3), we wish to define the same process 
decimated into even and odd terms giving 







C^z' 1 ) C 0 (z' 1 ) 




' y 0 m ' 




D(z~ l ) D(z ' 1 ) 




~e e [l] 


r— i 
1 

CN 

o 




z' 1 C 0 (z~') C e (z~ l ) 




_e„m_ 






1 

i 

<N 

1 





where 



e e [l] = e[2l], 


(5.9) 


o 

II 

l 

H— 1 


(5.10) 


D(z 2 ) = A(z' l )A(-z '), and 


(5.11) 


) + z~ l C 0 (z~ 2 ) = B(z~ ] )A(-z~ : ) , 


(5.12) 



where 



C e are the even coefficients of the expression B(z ')A(—z ') and 
C 0 are the odd coefficients of the expression B(z -1 )A(-.<: _1 ) . 

By breaking up the original signal y 0 [«] into even and odd equations, we get an 

identically equivalent system of two equations driven by two independent white noise 
processes at one-half the sampling rate of the original signal. If in effect we can 
construct the sequence y 0 [/i] from two independent white noise processes, then by the 
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properties of the ARMA model described in Chapter 2, the inverse processes will take the 
sequence y 0 [n] and create two independent white noise sequences. 

The proof for the even/odd decomposition of a signal shown in equation (5.8) 
starts with the ARMA filter determined from the original signal y 0 [n] such that 



H(z") 



B(z -') 



(5.13) 



In order to apply the Noble Identity later on in the proof, all transfer functions will 
be described in terms of z’ 1 instead of z . We want to represent the ratio of polynomials 
in (5.13) by separating it into even and odd powers of z' 1 giving us 

H(z~') = H e (z~ 2 ) + z~'H 0 ( Z - 2 ). (5.14) 

To determine what H e (z _1 ) and H 0 (z~') are, we assume that both transfer 

functions have the same denominator polynomial as a function of z" : • We accomplish 
this by multiplying the numerator and denominator by A(-z ~' ) , i.e. 



H(z' l ) = 



A(z -1 )- A(-z"') 



C(z~‘) 
D(z~ 2 ) ’ 



(5.15) 



Recall that a polynomial can be decomposed into even and odd parts as in 
equation (5.14). Performing this decomposition on the numerator C(z"') we get the 
polyphase decomposition of H(z~ l ) as 



H(z-') = H e ( Z - 2 )+z-'H 0 (z- 2 ) = 



C,(z' 2 ) 

D(z~ 2 ) 



C 0 (T~ 2 ) 
D(z~ 2 ) ' 



(5.16) 
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The numerator of equation (5.16) can be thought of as separating all the even and 
odd powers of z _1 from the multiplication of polynomials B(z ~ ') and A(-z -1 ). To 
show the denominator is a polynomial of powers of z~ 2 , however, takes some 
explanation. The denominator is the multiplication between polynomials A(z -1 ) and 
A(— z _1 ) . We choose to factor the polynomial A(z _1 ) in terms of its poles p t as 

ACz-VfK 1 -^- 1 ), (5.17) 

;=i 

where the total number of poles is also the filter length P . 

In the same way, A{—z ~ x ) can be factored as 

A(-z-') = f\(l +Pl z-'). (5.18) 



Substituting equations (5.17) and (5.18) into the denominator of the polyphase 
decomposition equation (5.16) we get 




(5.19) 



Again, the single series product reflects the fact that the same index is used for 
both expressions. You can see that the cross terms of the expression inside the series 
product are zero so that the resulting expression yields 

A(z'')- A(-z _1 ) = ]^[(1- /7, 2 z" 2 ) = Z)(z" 2 ). (5.20) 

;=i 
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From the ARMA decomposition equation (5.8) we see that a signal can be 
synthesized from white noise processed through even and odd parts of the ARMA model 
as in Figure V-3. 




Figure V-3 Block Diagram of ARMA Decomposition Equation 



Now consider each output in Figure V-3 separately. Starting with the even terms 
y 0 [2l] of the original signal we have the block diagram as shown in Figure V-4. Note 

that the downsampling operation is taken to each branch on the other side of the 
summation. At this point we use the Noble Identity as in Figure V-5 and we get the even 
and odd parts of the ARMA model in terms of z~ l shown in Figure V-6 [16]. 
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Figure V-4 Block Diagram for Even Terms of Original Signal 
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Figure V-5 The Noble Identity 




Figure V-6 Block Diagram for Even Terms After Applying the Noble Identity 
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From Figure V-6 we see that the white noise process e[n ] is split up into its even 
and odd terms as shown in figure V-7. 




y 0 m 



Figure V-7 Block Diagram to Determine y 0 [ 21] From a White Noise Process 



The equation to find the even terms of signal y 0 [»] is then 

>’ 0 [2/] = H e (z~ l )e[2l] + H a (z~')e[2l -1]. (5.21) 

Next we complete the proof by considering the odd terms of the ARMA 
decomposition equation (5.8) as in Figure V-8. Following the same concept used to 
compute the even terms of y 0 [n], we apply the Noble Identity again leading to the result 

shown in Figure V-9. 
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Figure V-8 Block Diagram for Odd Terms of Original Signal 




Figure V-9 Block Diagram for Odd Terms After Applying the Noble Identity 




Figure V-10 Block Diagram to Determine y 0 [2l - 1] from a White Noise Process 
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Therefore, Figure V-10 shows that the expression for the odd terms of signal 
v 0 [n] is given by 

y 0 [2l - 1 ] = H e (*-' )e[2l - 1] + z' 1 H 0 (z H )e[2l] . (5.22) 

Now consider the even and odd terms of the white noise process. Point by point, 
there is no correlation between the two processes and therefore, it can be shown that a set 
of odd and even terms of a single white noise process is completely independent of each 
other. To reflect this uncorrelatedness, we re-define the odd and even parts of a single 
white noise process as two completely independent white noise processes as 

e e [l] = e[2l] , and (5.23) 

e 0 [l] = e[2l-l], (5.24) 

It follows then that the complete model of the signal y 0 [/z] is the combination of 

its odd and even parts as defined in equations (5.21) and (5.22) and Figures V-7 and V- 
10. Combining these results along with the newly defined independent white noise 
processes e e [l ] and e Q [l], we get the system of equations 

y 0 [2/] = H e ( Z -' )e e [/] + H 0 (z’ 1 )e„ [/] (5.25) 

y 0 [2/ -1] = H e (z- l )e 0 [l]+z-'H 0 (z-')e e [l] , (5.26) 

and a block diagram described in Figure V-ll. Finally, by combining equations 
(5.25) and (5.26) into matrix form we get 

' y 0 [ m 

_y 0 C2l-l] 



‘ H e (z-') 




\ll] 


_z-'H 0 (z -') 


H e (z-')_ 


_e 0 [l]_ 



(5.27) 
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which is the polyphase ARMA of equation (5.8). 




Figure V-ll Block Diagram to Model a Signal from Two Independent White Noise 

Processes 
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We have shown that it is possible to model the even terms of a signal separately 
from the odd using two independent white noise processes. Based on this fact, along 
with the ARMA model given from the original signal v 0 [«], we wish to determine the 

corresponding ARMA model for y,[/] . Recall y^/] is the lowpass approximation of the 
original signal y 0 [rc] as described in the multiresolution operation of Figure V-2. 

As it has been said earlier, the multiresolution ARMA filter bank is dependent on 

the filter F(z _1 ). So from the averaging filter F(z _1 ) = -^-(1 + z~ l ), vj/] becomes 



>’,['] = ^y 0 [2/]+|y 0 [2/-u 



'i r 


y 0 [2/] ‘ 


2 2 


_y 0 [2/-i] 



(5.28) 



Substituting equation (5.27) into (5.28) we get 



>',[/] = 



‘i r 




H 0 (z- ] )' 


~e t [l] 


2 2 









(5.29) 



Carrying out the multiplication and substituting H e (z ')and H 0 (z~') in terms of 

c,(z"> c„(r‘) 



, and 
D(z - 1 ) D(z ) 



described in equation (5.16), we get 



ydn=- 



CAz^) . Cjz") 
D(z~ l ) ^ D(z~ l ) 






CM ) , c a (z ") 

D(z~ l ) D(z“) 



e o m 
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Making one final substitution the final expression becomes 




(5.30) 



where 



FM~ l ) = \[cAz- l ) + z- l C 0 (z~ 1 )\ and 

F 0 (z~') = \[c e (z~ l ) + C o (z-')}. 



(5.31) 



(5.32) 



Now consider equation (5.30) as a difference equation. We want to determine 
[/] as it satisfies the equation 



where 

d p = coefficients to the polynomial expression D(z~ l ) assuming d 0 = 1, 
f n e = coefficients to the polynomial expression F e {z~ x ), 
f° = coefficients to the polynomial expression F 0 (z”‘), 

P = length of original Autoregressive (AR) model A(z -1 ), 

Q = length of original Moving Average (MA) model 5(z _I ) , 

M =2P-l, and 
N = P + Q-l. 






N 



M 






(5.33) 



n=0 
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In order to solve this difference equation, we convert equation (5.30) from a sum 
of two transfer functions into a sum of two state space models. Recall the Type I 
controllable canonical state space form for the Direct Form realization of the transfer 
function [13] 



is 






b n z A +b.z r ' + ■ 



N - 1 



■K 



z N +a,z N ~' +••■ + «„ 



(5.34) 



1 

jc: 

sT 

+ 
H— ^ 

1 




0 


i 


v 2 [n + l] 




0 


0 


* 




0 


0 


v N [» + !]_ 




r a s 


~~ a N - 1 



0 



0 

0 

1 

- a, 





v,[n] 




'o' 




v 2 [n] 




0 


* 




+ 










1 



x(n) 



(5.35) 



Therefore, equation (5.35) is of the form 

V[n + 1] = AT[H] + B4«]> 

and 



y[*0 — [(^at b 0 a N ),(b N _] ^o^iV-i ), '”, (^i Vi)]' 



v, [«] 
v 2 [»] 






which leads to 



y[n] = C-V[n] + Dx[n]. 



(5.36) 
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Since there are two transfer functions in equation (5.30), we apply the state space 
form twice yielding 



A-J/ + l] = A, x[Z] + £,<?[/] 


(5.37) 


y u [/] = C r *,[/J+£>,^] 


(5.38) 


aJ/ + 1] = A ( ,a[/] + 5^[Z) 


(5.39) 


y,. 0 m = c 0 x 0 m+D 0 e[/] 


(5.40) 


y\U] = y Ue U) + y Uo [i], 


(5.41) 



where 



e[l] = 



e,m 

e 0 U] 



A e =\ = 



0 1 
0 0 
0 0 



a N a N - 1 



B,=B 0 = 



0 

0 

1 

~ a i 



C'=[{f' N -/‘o^),(/Vi ),•••, (/', -Tod,)] 

c a = [(/% -f°od N ),(/Vi -/V*., ),•••, (/°i -rod,)] 

D e = / V and 

D a =f°o. 
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A t = A a because these matrices are based on the same denominator term D(z 1 ) . 

From the state space model we can now express equation (5.30) in terms of matrix 
variables as 



Note that D e and D a are scalars. There are clear similarities between the two 
terms in equation (5.43), which is advantageous because we wish to combine the two 
expressions into one state space model. Since A e and A 0 are the same matrix, the term 
to determine the eigenvalues is equivalent. We combine both terms of equation (5.43) 
resulting in a single state space model for v,|7] , giving 




(5.42) 



Taking the transpose of expression (5.42) leads to 



* M = [bJ ( zI-A/y 1 cj + D e ]e e [l) + (bJ o zl-A/y 1 C 0 T +D 0 )e o [l}. (5.43) 



x[l + 1] = Ax[l] + Be[l] 

y l V] = Cx[l] + De[l], 



(5.44) 



(5.45) 



where the matrices in the new state space model are defined as 
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Looking at the new single set of state space equations (5.44) and (5.45), they 
describe a general stochastic state space model 



An innovation model can be determined from equations (5.46) and (5.47) using a 
Kalman Filter as 



z[k] is the observation at time k, 

a[A'] is the prediction of state x[k] given all observations up to time k - 1 , and 
z[k] is the innovation, which is independent on all past observations. 

The Kalman gain matrix K and the covariance matrix P are computed from the 
Riccati equation where 



x[k + l] = Ax[k] + v[k] 
z[k] = Cx[k] + \£k]. 



(5.46) 



(5.47) 



*[fc + l] = Ax[k]+ K-z[k], 

z[k]=z[k]-Cx[k]> 



(5.49) 



(5.50) 



where 



P = £{.v[£]3c r [A']} 
with x[k] = x[k]-x[k]. 



Finally, the covariance of the innovation is given by 
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E\nkf\=CPC' + R 



where R is the variance of the white noise process vv[&] . 



For the particular case of equations (5.44) and (5.45), we obtain 

x[l + l ] = Ax[l] + K-e 1 [l], (5.51) 

>\U] = Cx[l] + e ] [l). (5.52) 



In order to solve the Riccati equation, covariance matrices Q , R and S must be 
determined. Because we assume zero-mean white noise, the covariance matrices are 
identically equal to the correlation matrices. Applying the general state space model 
equations (5.46) and (5.47) with the specific model equations (5.44) and (5.45), we define 
Q as 

Q = e\v[]c] • v T [*]}= E{Be[l]e T [l)B T }. 



Performing the product of the uncorrelated 2 x N white noise set, we note that 



eU]e T U] = 



e e [l] 

e 0 U] 



km e 0 [l)} = 



CT‘ 0 
0 cr 2 



= o-I. 



Now substituting back into the above equation and noticing that all the 
expressions in the expectation are deterministic we get 



Q = o 2 BB T . 



(5.53) 
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Performing similar operations to R and S we end up with 



R = £■{ w[/]_w r [/]}= a 2 DD T , 



(5.54) 



and 



S = f{v[/] w T [I)}=a 2 BD T . 



(5.55) 



Note here that the variance <j 2 is the same in equations (5.53) through (5.55) 
because the white noise processes that comprise e[l ] are independent, yet have the same 
variance. Recall this variance is determined from the initial ARMA model. 

For the final step in determining the multiresolution innovation process, consider 
the measurement equation (5.52) where we have determined all values in order to 
determine the innovation e x [l ] as 



It would be possible to implement this multiresolution innovation process in terms 
of a Kalman Filtering recursion, but notice the state space model representation of 
equations (5.51) and (5.52). This set of equations is also in Type I controllable canonical 
state space form for the Direct Form realization of a transfer function [13]. This means 
that a single transfer function can be determined from the state space model, giving 



e l U) = y\U]-Cx[l) 



where it can be shown that the variance of e,[/] is 




(5.56) 



y l [l) = [C(zI-A)- l K + l]e l [l]. 



(5.57) 
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Finally, it can be shown that if an ARMA model is minimum phase, then all its 
corresponding multiresolution ARMA filters are minimum phase as well. From chapter 
2, recall that a minimum phase filter will have a causal stable inverse. From this fact, we 
see that the inverse to the multiresolution ARMA model yields our goal of obtaining a 
multiresolution innovation filter bank. 

In a similar way, we can determine a model for the highpass component, as in 
Figure V-l. Consider now the ARMA filter coefficients from the highpass filter channel 



where G(z ') = 



2 



-l 

2 



We see that the derivation is quite similar. In fact, the only 



difference is the section where we consider the effect ofF(z _1 ) on the Direct Form 
realization of the transfer function in equation (5.42). Going back to equation (5.29) we 
substitute the highpass filter G(z _1 ), giving 






'1 -f 


‘ H'{z~ x ) 




\[l] 


.2 2 _ 


_Z- l H 0 (z- 1 ) 


H e (z~ l )_ 





(5.58) 



Carrying out the multiplication, substituting in equation (5.16), and grouping 
terms yields the expression 



y.P] 



G'{z~ x ) 

D(z~ l ) 



e e [l] + 



G 0 (z~ l ) 

D{Z' X ) 



e 0 U], 



where 

G t {z") = \[cM~ l )-z x C 0 {z- x )] 



(5.59) 



(5.60) 



and 
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(5.61) 



G.<z' i > = tIc.(*' , )-C.U'‘)]- 



From this point, we can substitute G e {z ') and G 0 {z ') into equation (5.42), 

which is the state space form for the Direct Form realization of the transfer function in 
equation (5.59). 



B. METHOD 2 SUMMARY 

Method 1 was a method to look at a whitening filter before wavelet 
decomposition. By contrast, Method 2 investigates a method to perform the whitening 
filter after the multiresolution process. The algorithm can be divided into two processes. 
The first one is the standard wavelet decomposition of the original signal similar to the 
DWT on the whitened signal for Method 1. The second process in the algorithm for 
method 2 deals with determining the innovation filter bank. A set of ARMA filter 
coefficients is necessary for each wavelet decomposition channel, however all that is 
needed a priori is the ARMA model coefficients of the original signal. That is to say, we 
determine the ARMA filter coefficients based on the original signal and from that model 
we can determine all follow-on multiresolution channels. From this description, we can 
see how the Kalman filtering procedure described in this chapter is applied and a short 
description of the steps for the entire process is shown in Table V-l. 

All the white noise channels are uncorrelated with each other resulting in an 
ability to observe the innovation process with varying degrees of highpass and lowpass 
filtering. By way of analogy, one of the useful properties of the wavelet decomposition is 
denoising. Since most useful signal information occurs in lower frequencies for 
underwater acoustic applications, the wavelet decomposition can leave out some highpass 
information (i.e. details) upon synthesis thereby removing noise from the signal. 
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In our particular application, a similar method of running a signal through a 
constant-Q filter takes place, only now we take each channel and process it through a 
whitening filter. In effect, each channel has less and less high frequency content. It 
would then stand to reason that this multiresolution whitening filter bank might perform 
two desirable functions. First, the filter bank might be able to distinguish low frequency 
narrowband signals well due to its constant narrowing low pass filter, providing high 
resolution in low frequencies. Second, since the channels are uncorrelated, indication of 
a transient signal might be evident in one channel and not in others, depending on its 
frequency content. 
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STEP 


DESCRIPTION 


EQUATION 


MATLAB FUNCTION 


1 


Determine ARMA coefficients of original 
signal and variance cr 2 


2.3 


[b, a] = cristi _ ar( y ) 


2 


Decompose a filter transfer function into 
even and odd parts for Low Pass Filter 
F(z”) 


5.15,5.16 


[C t ,C 0 ,D]= polyd(num,den ) 


3 


Determine F e (z~' ) and F 0 (z _1 ) 


5.31,5.32 


f 


l^n+l > ^n+1 > a n+\ ’ a n + 1 J — 

nodel2 (b n , a n ) 


4 


Determine the Type I Controllable 
Canonical Form for H e (z~ l ) and// 0 (z"‘) 


5.35, 5.36 


[a„b„c„d,)= 

t/2ss(F„D) 


5 


Combine even and odd processes to one 
state space model 


5.44, 5.45 


r 


Pn+\ j i a n + 1 > a n+l J 

nodel2(b n ,a„ ) 


6 


Determine Q , R and S 


5.53-5.55 


Q = cf 2 BB t , R = <7 2 DD t 
S = (7 2 BD t 


7 


Determine Kalman Gain and State 
Covariance Matrix from Riccati Equation 


[14] 


[K,P] = 

mydlqe(A,G,C,Q,R,S) 


8 


Determine new variance cr 2 


5.56 


a 2 =CPC T +R 

i 


9 


Determine next multiresolution filter 
transfer function 


5.57 


\tuun,den\ — 
ssltf (A, AT,C,1) 


10 


Repeat steps 2 - 12 for HPF G(z~ l ) 


5.60, 5.61 


Included in Step 5 


11 


Perform 2 - 13 to the desired number of 
stages 




me th _0 102 


12 


Run the signal through the multiresolution 
whitening filter bank as in Figure V-l 


Figure V-l 


meth_0102 



Table V-l Process to Determine the Multiresolulion Whitening Filter Bank 
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VI. RESULTS 



Now that we have set a framework for two methods of multiresolution filtering, it 
is necessary to compare the techniques against the benchmark signal described in Chapter 
1. Method 1 can be described as pre-whitening a signal before multiresolution analysis, 
which is the Discrete Haar Wavelet decomposition. By contrast, a signal processed 
through Method 2 uses the same two processes in opposite order. Again, the benchmark 
signal is a fan blowing, which is representative of the colored noise environment of the 
heavily - trafficked coastal waters because it contains more of the lower frequencies due 
to a rotating blade through a fluid. 

The test signal used in this section is the first three seconds of the benchmark 
signal unless otherwise specified. The assumption is there are no other transients in this 
regime except for the synthetic signal that starts after the first second. A wide range of 
synthetic transients is tested from narrowband single frequency sinusoids to broadband 
white noise. The object of this section is to compare established detection techniques 
with Methods 1 and 2 using a wide range of transients. 

Since the STFT is the current standard for underwater acoustic analysis, all results 
will begin with a spectrogram followed by the whitening filter process. In this way, the 
viewer can receive an intuitive “feel” for comparing methods. Introducing the new 
detection methods, as we shall see, results in a higher confidence in positively identifying 
the transient. 

Finally, for ease of comparison, Method 1 and 2 results will be presented 
concurrently at each of the three transient tests. Because of the multiresolution process, 
each method has multiple channels to display. Since not all channels will indicate the 
presence of the signal, all figures left out are assumed to be clutter. For convenience, a 
small diagram of each process will precede every figure with its respective filter path 
highlighted. Knowing which channel is being displayed is important when analyzing the 
effect of the filtering process on the transient signal. 
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A. 



BENCHMARK TEST 



The method comparison will test three different transient types. The first two are 
narrowband transients at opposite sides of the frequency spectrum while the last test will 
be the more general broadband type. We will start with examining the narrowband 
transient signals and make adjustments to Methods 1 and 2 as necessary. 

1. Narrowband Transient: High Frequency Sinusoid 

Note that the maximum testable frequency is half the sampling frequency (just 
over 5 kHz) as the benchmark signal is sampled at 11.025 kHz. Here, we arbitrarily 
choose a 5 kHz sinusoidal transient with a duration of .1 second. The spectrogram in 
Figure VI-1 shows the transient is emerging from the background but is barely 
recognizable audibly because the amplitude is small. From this signal representation it 
can be seen that since the benchmark contains little frequency in the region above 2.5 
kHz, a transient signal above this frequency can be easily recognizable, assuming the 
amplitude is high enough. 

Figure VI-2 shows a simple diagram of Method 1 in order to keep track of the 
multiple outputs in the analysis tree followed by Figure VI-3, which shows the innovation 
process on the benchmark signal. Again, in this representation the transient is noticeable 
between the vertical dashed lines. Notice, though, the mean and variance of the noise 
floor compared to Figure VI-5, which is the first stage of the multiresolution process. In 
this procedure, the whitened signal is now processed through a lowpass filter and a 
highpass filter, which are the Haar Wavelet coefficients. 

Looking at the results of the first stage of the Haar Wavelet decomposition shown 
in Figure VI-5, the signal through the lowpass filter provides no useful information, but 
the highpass filter shows the transient very clear and w'ith less variance on the noise floor. 

As a minor digression, it is important to recall that the signal passed through the 
Haar Wavelet filters will result in the Haar Wavelet coefficients, which satisfy Parseval’s 
theorem similar to the properties of the coefficients determined from the Fourier 
Transform. In this application, the actual numerical power given for each scale is not 
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important, only their relative values and so all the signals can be normalized. For this 
application, the fundamental properties of the Discrete Wavelet Transform are only 
important in terms of defining the filtering process. To effectively analyze the results, we 
simply view the multiresolution filtering process in terms of how the filters effect the 
signal in the frequency domain. What we have shown in Figures VI-3 and VI-5 is a 
denoising process, which is one of the many uses of wavelet transforms. 
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Figure VI-1 Spectrogram of Transient 1: 5 kHz Sinusoid 
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Figure VI-2 Method 1 Analysis Tree (\v[n]) 




Figure VI-3 Transient 1: Whitening Filter 
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Figure VI-4 Method I Wavelet Analysis Tree (A, and £>,) 




Figure VI-5 Transient 1: Method 1 , Stage 1 (A, and £>,) 
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For the high frequency sinusoid, the transient will not appear in any follow-on 
stages of the lowpass filtering processes (approximations). At this point, we are left with 
a successfully detected signal in the high pass channel of the first stage. For this 
scenario, we can see the transient above the noise floor, but recall this signal can also be 
detected from the spectrogram. Now, suppose we can reduce the amplitude of the 
transient signal to the point where it cannot be easily identified in either the spectrogram 
or the whitening filter. Can the denoising process of Method 1 provide better results? 
The answer is shown in the next set of figures. 

The key to successful denoising on the high frequency end of the spectrum is to 
filter the details (high pass channel) in the same manner as the approximations (low pass 
channel). Knowing that we can remove the signal from the noise on the low pass 
channel, we apply the same wavelet decomposition at every output. This description is 
known as the standard Wavelet Packet framework, which has been applied to a wide 
range of applications such as denoising and data compression. Wavelet packet 
processing allows for a more complex and flexible analysis because both the details 
(highpass channel) and the approximations (lowpass channel) are decomposed, as shown 
in Figure VI-7. In terms of wavelet analysis, the wavelet packet framework allows for a 
wider range of bases from which you can choose the best representation [15]. 

For our application, the flexibility of the wavelet packet allows for greater 
narrowband detection sensitivity in both high and low frequencies. Take the high 
frequency transient for example. Let’s reduce the amplitude of the 5 kHz sinusoid from 
.005 to .003. To quantify this reduction, audibly detecting this transient signal is like 
straining your ears to listen for the softest high frequency sounds during an audiogram 
test. On top of that, we still have the noise of the fan blowing. Figure VI-6 shows the 
spectrogram of the benchmark signal with the barely noticeable transient at 5 kHz, which 
could be classified as noise. Now we add the high pass channel decomposition and find 
that several channels have successfully reduced the noise floor so that the signal can be 
detected. Because of the branching effect of the wavelet packet, a three-stage framework 
will yield 15 plots. In order to reduce the number of plots, only the graphs that detect the 
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signal are plotted. Figure VI-8 shows the whitening filter representation of the 
benchmark signal. Notice that the signal is completely imbedded within the noise such 
that it cannot be identified. Performing the first stage of the wavelet decomposition as 
before yields a slight indication of the transient (Figure VI- 10), but now we perform the 
filtering process on the details and notice the signal clearly emerges from the noise floor 
(Figure VI-12). Performing the filtering process one more time yields an even better 
result (Figure VI- 14). Realistically, the transient might have been positively identified on 
the spectrogram. However with the wavelet denoising process, the added signal detection 
capabilities might add to the measure of confidence. 

At first glance, it would appear that the 5 kHz signal should be detected in the bin 
to the far right marked DDD 3 because the bandpass filter in this region contains the 
frequency of the transient. Note in Figures VI-11 to VI-14 that the signal appears in the 
bandpass filter region AD 2 and AAD 3 . The reason for this could be due to aliasing that 

comes from the downsampling process is subject for further research. Also, since the 
highpass and lowpass filters of the Haar wavelet do not have a sharp dropoff at the cutoff 
frequency, it is possible to have part of a high frequency transient signal to appear in the 
low frequency channel and vice versa. Both these issues are subject for further research. 
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Figure VI-6 Spectrogram of Modified Transient 1: Reduced Amplitude 
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Figure VI-7 Wavelet Packet Framework (w[n]) 




Figure VI-8 Modified Transient 1: Whitening Filter 
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Figure VI-9 Wavelet Packet Framework (A, and £>,) 
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Figure VI-10 Modified Transient 1: Stage 1, Method 1 (A, and D,) 
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Figure VI-11 Wavelet Packet Framework ( AD, and DD, ) 
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Figure VI-12 Modified Transient 1: Method 1, Stage 2 (AD, and DD,) 
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Figure VI-13 Wavelet Packet Framework ( AAD 3 and DAD 3 ) 
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Figure VI-14 Modified Transient 1: Method 1, Stage 3 ( AAD - and DAD 3 ) 
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Turning now to Method 2, recall that the input signal is whitened after the 
multiresolution Wavelet decomposition as described in Figure V-l. Throughout the 
derivation, we used a slightly modified filter bank from the Haar Wavelet coefficients as 
described in equation (5.7). As mentioned before, the only difference is a constant factor 
and since the results are normalized, the shape of the signal is exactly the same. For the 
sake of consistency, the actual Haar Wavelet filter coefficients were used. 

Unfortunately, Method 2 did not detect the narrowband high frequency transient 
signal at amplitude .005. Looking at the lowpass filter, notice that it averages two 

adjacent data points and weights them by a factor of 4l . For the case of the 5kHz 
signal, perhaps the transient cannot be seen in the lowpass channel because of its high 
frequency content and cannot be seen in the highpass channel because the transient is 
averaged in with the noise. Nonetheless, to show that the method is still capable of 
detection, it can be shown experimentally that if the amplitude of the transient signal is 
increased to .01, the transient signal emerges from the noise. 

Although the multiresolution innovation process did not perform as well as the 
wavelet packet denoising process at 5 kHz, it is still a viable detection scheme. Indeed, 
we shall see other examples where Method 2 provides the best transient detection 
capabilities. 
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2. Narrowband Transient: Low Frequency Sinusoid 

At high frequencies, the spectrogram works quite well in detecting narrowband 
transients because the stationary colored noise model has very little high frequency 
content. Therefore any added frequency content is easily recognizable. Now consider a 
transient signal buried in the low frequency band between 0 and 300 Hz. Since the 
colored noise signal has its frequency content in this range, a signal cannot be 
distinguished from the background noise. Since both Methods 1 and 2 are based on a 
signal’s statistics in the time domain and not on the Fourier Transform in the frequency 
domain, they should detect low frequency transients equally as well. At each stage of the 
wavelet decomposition, the low frequency channel continues to be passed through a 
lowpass filter, effectively cutting the frequency content of the remaining signal in half. 
Continuing this process provides infinite resolution in the low frequency channels, which 
again is one of the nice properties of the Discrete Wavelet Transform. For this reason. 
Method 2 should perform at least as well as Method 1 in the high frequency narrowband 
case. 

We choose an arbitrary sinusoidal signal at 50 Hz and conduct a similar analysis. 
Figure VI- 15 shows the spectrogram of the benchmark signal along with the imbedded 50 
Hz transient signal. It is hard to tell whether or not the transient can be detected in this 
region because it could easily be identified as a spurious point. The transient clearly 
cannot be detected with just the whitening filter (Figure VI- 17) and not even audibly. 
Starting with Method 1, we see the signal emerging on the low pass channel in Figure VI- 
19. Because the frequency is so low, the signal becomes more pronounced with each 
successive lowpass filter as shown in Figures VI-20 through 23. 
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Figure VI-15 Spectrogram of Transient 2: 50 Hz Sinusoid 
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Figure VI-16 Wavelet Packet Framework (w[n]) 
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Figure VI-18 Wavelet Packet Framework (A, and £>,) 
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Figure VI-20 Wavelet Packet Framework ( AA 2 and DA 2 ) 
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Figure VI-21 Transient 2: Method 1, Stage 2 ( AA 2 and DA 2 ) 
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Figure VI-22 Wavelet Packet Framework ( AAA 3 and DAA 3 ) 
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Figure VI-23 Transient 2: Method 1, Stage 3 (AA4 3 and DAA 3 ) 
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Method 2 results for the low frequency signal are just as revealing. Unlike the 
Wavelet Transform on the white noise process, the colored noise signal does not have 
equal frequency content at all frequencies (except for the transient) and therefore the 
narrowband transient will not fall on one side of the highpass/lowpass filtering process. 
Indeed, the signal may show up in as many as all of the channels as is the case here, 
which in and of itself might provide added redundancy since the channels are 
independent of each other. For the lowpass-filtering channel, the benchmark signal is 
downsampled and averaged. As it averages adjacent points, the process essentially 
reduces the noise level while maintaining a strong indication of the presence of a signal. 
For the highpass filter, the difference between adjacent points is computed and the 
combination of the two processes form a basis set for the original signal. Again, this 
highpass process might prove valuable since it is independent of the lowpass process. 
Figure VI-25 shows the signal clearly emerging in both channels. With each successive 
multiresolution innovation process, the noise floor is reduced and the signal becomes 
more pronounced (Figures VI-26 to 29). 
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Figure VI-24 Multiresolution Innovation (A, and D x ) 




Figure VI-25 Transient 2: Method 2, Stage 1 (A, and D x ) 
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Figure VI-26 Multiresolution Innovation Stage 2 ( AA 2 and DA 2 ) 
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Figure VI-27 Transient 2: Method 2, Stage 2 ( AA 2 and DA 2 ) 
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Figure VI-28 Multiresolution Innovation Stage 3 ( AAA 3 and ZMA 3 ) 
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Figure VI-29 Transient 2: Method 2, Stage 3 ( AAA^ and DAA 3 ) 
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3. 



Broadband Transient: Elvis 



We have studied the effects of Methods 1 and 2 on narrowband transients. In the 
case of underwater acoustics, it is quite possible that the signal we wish to detect is a 
single frequency like an active sonar pulse or a unique single frequency tonal. On the 
other hand, a broadband signal such as the splash of a mine hitting the water or the clang 
of a metallic object is also possible. 

From the narrowband transient, we now consider a more general case where the 
expected signal has a wide frequency band. For this test, the frequency content must be 
completely arbitrary, so for no apparent scientific reason the arbitrary colored noise 
signal will be Elvis singing that great Christmas classic, “Blue Christmas’’. Besides, it is 
just as likely to find Elvis in the deep blue ocean as it is to find him anywhere else. 

In order to see the frequency content of the narrowband transient signal, a 
spectrogram of the transient alone is provided in Figure VI-30(a). Notice that Elvis can 
really hit those low notes. The strength of the frequency content is large in the range 
between 0 and 1 kHz, then tapers off as it approaches 5 kHz. Adding Elvis to the colored 
noise model as in Figure VI-30(b), notice that the signal can barely be seen, and that only 
because we have other spectrograms of the same data with which to compare. Without 
the redundant benchmark data, “Blue Christmas” is reduced to mere noise. We perform 
the innovation process on the benchmark signal in Figure VI-32 and, no doubt about it, 
Elvis is clearly alive. Performing the denoising on the innovation as in Method 1 
confirms this to be true, however the signal is less pronounced. Also notice that the 
signal shows up - in part - on both highpass and lowpass channels at the same stage. 
This indicates that the broad range of frequency content of Elvis’ golden voice is divided 
into highpass and lowpass regions. The problem with Method 1 here is that along with 
the reduction of noise with each filter stage comes the reduction of the signal power as 
well depending on its frequency content. In other words, the success of Method 1 is now 
dependent on frequency content. 

The signal had the same effect when processed through Method 2, only worse. 
Only the first stage provided any indication of Elvis’ location as in Figure VI-44. Here 
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again, the broadband signal could have been averaged in with the noise so that the 
innovation process could not distinguish Elvis from a fan. 
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Figure VI-30 Spectrogram of Elvis Without Benchmark Colored Noise 
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Figure VI-31 Spectrogram of Transient 3: Elvis 



102 




Figure VI-32 Wavelet Packet Framework (w[n]) 




Figure VI-33 Elvis Through a Whitening Filter 
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Figure VI-34 Wavelet Packet Framework (A, and £>,) 
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Figure VI-35 In Search of Elvis: Method 1, Stage 1 (A, and D,) 



104 




Figure VI-36 Wavelet Packet Framework (AA and DA 2 ) 



1 

0.8 
0.6 
0.4 
0.2 
0 

0 0.5 1 1.5 2 2.5 

1 

0.8 
0.6 
0.4 
0.2 
0 

0 0.5 1 1.5 2 2.5 

TIME (seconds) 



Figure VI-37 In Search of Elvis: Method 1, Stage 2 ( AA Z and DA ) 
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Figure VI-38 Wavelet Packet Framework ( ADA l and DDA 3 ) 



LOWPASS FILTER (APPROXIMATION) 




0 0.5 1 1.5 2 2.5 

TIME (seconds) 



Figure VI-39 In Search of Elvis: Method 1, Stage 3 (ADA 3 and DDA 2 ) 
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Figure VI-40 Wavelet Packet Framework ( AD 2 and DD 2 ) 
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Figure VI-41 In Search of Elvis: Method 1, Stage 3 ( AD Z and DD 2 ) 
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Figure VI-42 Wavelet Packet Framework (ADD 3 and DDD 3 ) 
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Figure VI-43 In Search of Elvis: Method 1, Stage 3 ( ADD 3 and DDD 3 ) 
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Figure VI-44 Multiresolution Innovation Stage 1 ( A l and D x ) 
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Figure VI-45 In Search of Elvis: Method 2, Stage 1( A x and D, ) 
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B. 



SUMMARY OF OBSERVATIONS 



We have presented a representative set of transient signals tested on two detection 
schemes. This section aims to provide a broad brushstroke of the effectiveness of the 
Methods as they are compared to the spectrogram and the whitening filter. To 
accomplish this comparison, a measure of successfulness must be defined. Since the 
location and duration of the transient signal is known, an estimate of the signal power and 
the noise power can be calculated as the Signal-to-Noise ratio. Since our signal model 
has additive noise, the signal power is simply the power above the noise threshold 
divided by the noise power where both power terms can be calculated from the whitened 
signal. For Methods 1 and 2, whose multiple channels will not always provide the best 
transient detection, only the channel with the maximum SNR will be taken. The 
assumption here is that there is a way of always determining the best channel within a 
given method scheme. 

For a signal whose power is twice the noise power, the SNR will be 3 dB. This 
value will be our threshold detection mark where any SNR above this value is assumed to 
have positively detected the transient. Any value below this mark is assumed to not be 
able to identify the transient even though a visual representation might indicate 
otherwise. 

As in all experiments, the SNR measure is only an approximation. In our case, 
the same exact transient signals within the colored noise model is used to compute the 
SNR and so the results among the different methods may be compared to one another. It 
is fair to assume, however, that the results from this experiment are indicative of many 
types of transient signals hidden within an arbitrary colored noise scheme. The first 
experiment deals with broadband signal detection followed by a test on narrowband 
signal detection. 



110 



1. Benchmark Comparison 

The first experiment begins with Elvis. Taking the signal, we increase the 
normalized volume from 0 to 1 and determine its SNR from the whitening filter alone, 
Method 1 and Method 2. Figure VI-45 shows the results. Here, it can be seen that 
Method 1 reaches the 3-dB threshold first, followed closely by Method 2 and the 
benchmark whitening filter. As Elvis gets louder, though. Method 1 reaches a higher 
maximum SNR value of 25 dB while the other two methods are just about identical as 
they reach a steady state value of about 20 dB. The results of the experiments show that 
the noise level does indeed reduce faster than the signal and so a better detection method 
was found even for arbitrary broadband signals. On the other hand, comparing any of 
these results with the threshold reveals that Method 1 might reach detection status before 
the other methods, but not by much. What’s more, the spectrogram might be able to pick 
up the signal before any of the three signals, but again there must be some frequency 
content outside of the frequencies defined by the colored noise. 

Another broadband test was conducted where the frequency content expanded the 
entire spectrum. A similar test was conducted where the broadband signal was white 
noise and the results are found on Figure VI-46. Here again, Method 1 broke out with the 
highest maximum SNR at 40 dB. Since all three methods reached the 3 dB threshold at 
about the same amplitude, all three methods perform about the same. Also, we are 
guaranteed to have the spectrogram detect the transient at some point close to the other 
methods and so any of these processes are equally effective. 
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Figure VI-46 SNR Comparison on Colored Noise Transient 




Figure VI-47 SNR Comparison on White Noise Transient 
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The final comparison is with the narrowband transient. For this test we vary both 
frequency and amplitude of a standard .1 second transient and process this signal through 
the whitening filter, Method 1 and Method 2. The resulting 3-dimensional image reveals 
some interesting peaks and valleys. Figure VI-47 shows the result from the whitening 
filter. Comparing this plot with the results from Method 1 (Figure VI-48) we see right 
away that it is the same basic shape, only the SNR has increased by about the same 
amount over the entire frequency-amplitude plane. Now recall that the threshold is 3 dB 
and notice how well the wavelet denoising process worked. In the same manner, Method 
2 produced very good results in the low frequency region. 

We wish to further quantify the results of the comparison test. Taking the 3- 
dimensional plot, we slice the graph along the frequency-amplitude plane at a constant 
level of 3 dB. Looking down on the slice, we can see clearly the regions where the SNR 
is above 3 dB, which is the detection region. Figure VI-49 shows the modest narrowband 
transient region for the whitening process. With the exception of a small strip, the lower 
frequency SNRs remain below the threshold, especially between 0 and 300 Hz where 
there is no SNR above 3dB. Now Comparing Method 1 to the whitening filter detection 
region (Figure VI-50) we notice right away the larger areas in which the SNR is above 
the threshold. To quantify, Table VI-1 shows that only about 30% of the entire 
frequency-amplitude plane is above the threshold. With the Haar Wavelet denoising 
method, the percent detection over the entire region increased to 74%. In the low 
frequency region between 0 and 2 kHz, the detection region increased from 29.7% to 
78.3%. The results from Method 2 as seen in Figure VI-51 show an equally strong 
performance in the low frequency regions and as mentioned earlier, it has limitations in 
the high frequency regime. 

These results seem very good compared to the whitening filter, now compare the 
methods against the spectrogram. The spectrogram can detect most signals above an 
amplitude of .1 and frequencies above 300 Hz. Recall that the whitening filter did not 
have any of its detection region below 300 Hz so with these two methods combined, low 
frequency transients would go unnoticed. Looking back at Figures VI-50, VI-51 and 
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Table VI- 1, we notice that the detection region above the 3 dB mark is now almost 50% 
of the total area in this region. 
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FILTERING 

TECHNIQUE 


PERCENT DETECTION IN FREQUENCY RANGE 


0 - 5 (kHz) 


0 - 2 (kHz) 


0 - 300 (Hz) 


WHITENING FILTER 


30.3 


30.3 


0.0 


METHOD 1: HAAR 
WAVELET DE- 
NOTING 


73.9 


79.6 


50.5 


METHOD 2: 

MULTIRESOLUTION 

INNOVATIONS 


66.1 


81.1 


51.8 



Table VI-1 Narrowband Transient Detection Capability Comparison Over a 

Threshold of 3 dB 
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Figure VI-48 Narrowband Transient SNR For Whitening Filter Process 
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Figure VI-49 Narrowband Transient SNR for Haar Wavelet Denoising 
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Figure VI-50 Narrowband Transient Detection Region For Whitening Filter 

(Threshold = 3 dB) 
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Figure VI-51 Narrowband Transient Detection Region For Haar Wavelet Denoising 

(Threshold = 3 dB) 
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Figure VI-52 Narrowband Transient Detection Region For Whitening Filter 

(Threshold = 3 dB) 
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Figure VI-53 Narrowband Transient Detection Region For Multiresolution 
Innovation Process (Threshold = 3 dB) 
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2. Strengths 

The most useful property of the wavelet denoising process and the multiresolution 
innovation process are their abilities to detect narrowband signals in the low frequency 
regime, particularly in between 0 and 300 Hz. This property is highly desirable since it 
has been shown that neither the spectrogram nor the whitening filter is effective in this 
frequency region. 

The Haar Wavelet denoising routine works well at high frequencies too because 
of its wavelet packet framework. Here, the high frequency details are filtered too, 
creating greater flexibility by isolating more frequencies. 

Since the multiple-channel outputs of both Method 1 and Method 2 are 
independent of each other, it is likely that more than one channel will indicate the 
presence of a transient signal. This property provides a measure of redundancy in 
positively identifying a transient signal. 

Both methods are simple to compute numerically; the DWT is computed by a 
simple filtering process as the number of computations is comparable to current methods 
that use the FFT. Since the size of the data set is reduced by a factor of 2 at each stage, 
the number of computations is reduced by the same amount. This property is desirable 
because we know that any follow-on stages will incur less and less computations. 

Both methods lend themselves to an efficient modular processing routine with a 
user-friendly GUI display. Take for example the wavelet packet framework. The user 
might designate the number of stages to explore and at any time, can add another stage 
without having to go back and compute the entire wavelet tree. The same is true for the 
multiresolution whitening process because, like the wavelet decomposition, every 
successive ARMA filter set is only based on the previous stage and does not need to be 
computed from the very first stage again. 



3. 



Weaknesses 



Both wavelet denoising and multiresolution innovation procedures attempt to 
reduce the noise variance by a series of filtering processes. For a signal that is isolated to 
a small range of frequencies, both methods will perform quite well since the filtering 
process has the ability to cut off a sizeable amount of extraneous frequency content. By 
contrast, a broadband signal will have part of its frequency content cut off at each 
filtering stage. The effect will be a reduction of signal power along with a reduction in 
noise power. Even with the signal power reduction, the results of the previous section 
show that Methods 1 and 2 were able to reduce the noise floor as the amplitude of the 
transient signal increased. Comparing these methods to the spectrogram, however, we 
see that the transient will appear as a vertical line at amplitude less than .05, which is 
about when the signal emerges from the noise in the multiresolution filtering process. 
The advantage of the spectrogram for broadband signals is its color representation, which 
essentially adds an extra dimension to detecting the transient. The only exception to this 
observation is when the transient contains all or part of the same frequencies as the 
background noise, in which case the signal is camouflaged in the spectrogram. 

Inherent in the wavelet decomposition process is a reduction of resolution in the 
time domain. As described in Chapter 3, to perform a DWT on a signal it is 
downsampled as part of the process. By this action, the signal in the time domain loses 
some detail at each stage. Figure VI-53 shows the entire 6 seconds of the benchmark. 
Recall the peaks located in the latter half of the benchmark signal are a set of 1 1 
broadband transients. Notice the cluster of pulses just before the fourth second in the 
graph of the highest sampling rate. Here, the three spikes are easily recognizable. Now 
dropping down through the other three graphs where the sampling frequency is reduced, 
notice that the three spikes are slowly combining to form a single transient. If the final 
stages of method 1 or 2 were the only channels to pick up this set of transients, there 
would be no way to differentiate between them. 
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Figure VI-54 Reduction of Signal Resolution Due to Downsampling 
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VII. CONCLUSIONS 



This thesis introduced new transient detection techniques designed to provide 
sonar operators an enhanced capability to detect enemy acoustic signals. The methods 
presented in this thesis are based on the Autoregressive-Moving Average and Discrete 
Wavelet Transform. The difference is the order in which these two operations are 
performed. Method 1 first whitens the input signal using an ARMA model before 
denoising the whitened signal with a wavelet decomposition operation. Method 2 first 
applies wavelet decomposition before applying a bank of whitening filters to each 
wavelet channel. 



A. USW LITTORAL WARFARE APPLICATIONS 

The experimental results reveal a significant increase in signal detection using 
either proposed method when the transient is narrowband. This ability to localize 
frequencies might be useful for applications where the frequency of a potential transient 
is known a priori, such as for example in active sonar operations where the return signal 
is known to be within a Doppler shift of the transmitted frequency. Wavelet denoising 
has the ability to narrow its filtering process such that a specific channel can be observed, 
rejecting all other frequency content. Therefore, this process has the potential to increase 
the effective range of the sonar without having to increase the signal power. 

Another application to narrowband signal detection might be for passive sonar 
operations where one attempts to uncover a known frequency within noise of similar 
frequency content. In low SNR, we have shown that the spectrogram does not provide 
useful capabilities in this regime and so we rely on time-domain methods.' Examples of 
such signals might be specific tonals or blade rates produced by enemy platforms hidden 
within the noise of background shipping. 

In spite of the fact that a computer can successfully detect a transient signal over 
some given threshold, the sonar operator is still a vital member in the detection phase of 
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USW. Indeed, it is not unusual for the initial detection of an enemy submarine to come 
from a source other than an automatic alert function. Compared to the other warfare 
areas, human intervention is more vital in the detection phase because the pace is 
generally slower. With this in mind, an operator must have multiple displays and sensors 
at his disposal. It is believed that the methods tested in this thesis will not only increase 
detection capabilities, but will also provide enhanced display features. 

The main idea behind multiple display features at an operator’s disposal is to 
correlate data. Consider a Graphical User Interface (GUI) system where an operator is 
shown a wavelet packet framework diagram as in Figure VI-7 where the boxes represent 
the wavelet decomposition channels. The display feature could also link with existing 
features in the sonar suite to further enhance transient detection. 

Another way to enhance operator displays would be to include various intensity 
levels that somehow correspond to the SNR. Examples include full color displays or 
alarms with a range of volumes relative to the strength of a transient signal. 
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B. 



SUGGESTIONS FOR FUTURE RESEARCH 



As the wavelet decomposition can theoretically be earned on to infinite detail, so 
can the ideas of improving technology. Although this section is certainly not all- 
inclusive, several suggestions are provided for future research based on the results from 
Methods 1 and 2. 

1. Method 1 Improvements 

The DWT denoising process of Method 1 uses the Haar Wavelet, which is the 
most basic wavelet basis function that uses only a first order polynomial for the highpass 
and lowpass filter coefficients. As one can see from the frequency response of the filter 
coefficients (Figure III-9), there is a significant overlap between the two filters, which 
results in a significant amount of high frequency information contained in the low 
frequency channel and vice versa. Further research can be done to study other wavelet 
basis functions that have a sharper frequency response. To see a preliminary example of 
denoising using a different wavelet basis function, refer to Appendix B. 

Because we know that the wavelet decomposition produces independent channels, 
another area for follow-on research is to define a way to correlate data at each channel, as 
described in Appendix B. 



2. Method 2 Improvements 

Recall the successful improvements the dyadic filtering had on the narrowband 
high frequency transients. It is believed that the same success will occur for Method 2 if 
the wavelet decomposition were carried out on the highpass channels (details) as well. 
This would require further investigation on defining the whitening filter coefficients for 
these channels. In effect, Method 2 would then perform the full wavelet packet 
decomposition on the input signal and then run it through the modified innovation filter 
bank. Like Method 1, the modified Method 2 would result in 15 channels for the three- 
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stage wavelet decomposition. Again, this structure lends itself to a useful GUI display in 
which access to a large number of signal data is relatively simple. 

In addition to investigating the wavelet packet structure for the DWT-Innovation 
process, the same argument can be made about possible improvements using another 
wavelet basis set such as the Daubechies family. This process is a bit more complicated 
than changing the wavelet set in Method 1 because the innovation filter bank of Method 2 
is dependent on the wavelet filter coefficients. 

3. General Improvements 

The benchmark signal was a colored noise model with synthetic signals 
representative of coastal water acoustics. The main thrust for future study should be to 
use actual acoustic data in this regime. In this way, a more realistic conclusion can be 
drawn for this particular application. 

Results showed that the denoising process performed well for both Methods 1 and 
2. Perhaps the two methods could be combined such that we take the whitened signals 
from Method 2 and pass them through the DWT denoising process of Method 1. 

Looking at the benchmark signal, note that the three-second window performed 
well for the .1 second transient signal. However, this window size might not work well 
for signals of different length. If we wish to detect a range of transient signal durations, 
say as long as 6 seconds and as short as .01 second, further research must be done to 
investigate the effects of different window sizes with varying transient lengths. 

Finally, the last suggestion for future research would be to investigate methods for 
applying a neural network scheme to classify patterns of signals among noise in order to 
detect transients at a lower SNR. If there is a way to determine features that positively 
identify transient signals, then for USW applications it can be another valid tool to make 
the detection process more robust. 
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APPENDIX A. 



HAAR WAVELET PROPERTIES 
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A. 



NECESSARY TIME DOMAIN CONDITIONS 



This appendix shows that the Haar wavelet system meets all 8 necessary 
conditions for a valid wavelet system [9]. Starting with the wavelet decomposition tree 
in Figure III-6 and the MRA equation, we wish to determine the Haar coefficients of the 
PR filter bank, the scaling function and the shape of the Mother Wavelet. Based on the 
set of eight necessary conditions, we will show that the Haar system of equations is 
indeed a valid wavelet system. 

Recall that all parameters that define a wavelet system can be determined from 
scaling function. If we can somehow develop a parameterizing scheme to define the 
scaling coefficients h(n ) , we could determine any wavelet system. To this end, we begin 
by choosing an arbitrary filter size N . The filter size for the scaling coefficients 
h(n) must be even and must satisfy both the linear constraint 



5>(n) = V 2 , (A.l) 



and the — bilinear constraints which leads to 
9 



SW J =1. (A. 2) 



It is important to note that the filter h(n) is an FIR filter, which also implies 
compact support. This wavelet property is highly desirable as it aids in the time- 
localizing ability of the DWT and it also reduces the number of computations. 

N 

It can be shown that there are — - 1 degrees of freedom (parameters) in choosing 

h(n) that will guarantee the existence of a valid wavelet system [9]. Consider the 
simplest case where N = 2. This filter length means that there are zero degrees of 
freedom and by equations (A.l) and (A. 2), we get 
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/j(0) + h{\) = yfl , 



(A. 3) 



and 



h 2 (0) + h 2 (l) = l. 



(A.4) 



Although the second equation is nonlinear, it can be seen by inspection that the 
unique solution to this set of equations is: 



These are the Haar scaling function coefficients, which were developed in 1910, 
separate from studies in modem wavelet analysis. In 1992, Ingrid Daubechies showed 
that Haar’s work on this set of orthonormal basis functions are in fact the simplest set in a 
larger family of wavelet systems. From equation (3.19), we can determine the wavelet 
coefficients as 



With just the coefficients to the Haar scaling and Wavelet functions, it can be 
shown that this wavelet system meets eight necessary conditions [9]. 




(A. 5) 




(A. 6) 
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1 . 



Theorem 1: 



If cp(t) £ if is a solution to the MRA equation (3.14), and if j* (p(t)dt -A 0 , then 

5>(k) = V 2. (A. 7) 

n 

For the Haar Wavelet the above equation leads to: 

/z(0) + /i(1)=4= + -7= = ^2- 

V2 V2 



2. Theorem 2: 

//~ #>(t) A an if solution to the MRA equation (3.14) with J <p(t)dt -t- 1, then 

5>(t-o=5>w =1 , (a. 8) 

/ / 

with <A>(n + 2/zfc) * 0 for some k, then 

YjK 2n) = y £ j h(2n + l), 

n n 

where (3.30) may have to be a distributional sum 
satisfied, then (3.30) is true. 

For the Haar Wavelet the above equation leads to: 

/z(0) + /z(2) = /z(l) + /z(3) 



(A. 9) 

. Conversely, if (3.31) is 
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3. 



Theorem 3: 



If (p(t)is an Lr n L 1 solution to the MRA equation (3.14) and if integer translates 
of <p(t) are orthogonal as defined by 



| (p(t)cp(t - k)dt = ES(k) = 



\E k=0 1 
[0 otherwise J ’ 



(A. 10) 



then 



^h(n)h(n-2k) = S k 



1 k = 0 1 
0 otherwise J 



(A. 11) 



For the Haar Wavelet the above equation leads to: 

h(Qi)h(-2k) + h(l)h(l - 2k) 

/K0)/z(0) + /z(l)/z(l) = ^ + | = l for k = 0 
h(0)-0 + h(l)-0 = 0 fork<0 and &>0. 



Corollary 1: 

Under the assumptions of Theorem 3, the norm of h(n) is automatically unity, 
i.e., 






(A. 12) 



For the Haar Wavelet, corollary 1 leads to: 



/r(0) + /r(l) = — + — = 1. 
2 2 



Corollary 2: 

Under the assumptions of theorem 3, 



2>(2n) = 2>(2» + l)=-I- 

n n V 2 



For the Haar Wavelet, corollary 2 leads to: 



h(0) + h(2) = h(l) + /:(3) = — =• . 

V2 



4. Theorem 4: 

If (p(t)has compact support on 0<t<N-\ and if (pit -k) are linearly 
independent, then h(n) also has compact support over §<n<N-\: 

h(n) = Ofor n <0 and n > N -l . (A. 13) 

For the Haar Wavelet the above equation leads to: 

Haar Wavelet has compact support as h(n) = 0 

for n < 0 and n > 1 . 
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B. 



NECESSARY FREQUENCY DOMAIN CONDITIONS 



The last four theorems are necessary conditions in the frequency domain. To 
begin the analysis in the frequency domain, we look at the frequency response of the 
scaling coefficients h(n) . The Haar Scaling filter coefficients in the z - domain form the 
equation 



h o(z~')=-j=(l+ Z~'). 



(A. 14) 



Substituting e JK for z we get the frequency response of the FIR filter 



H 0 (co) 




+ e' jw 



)■ 



(A. 15) 



Likewise, the frequency response of the Haar Wavelet filter coefficients is 




(A. 16) 



It is important to note here that the frequency response of the scaling coefficients 
h{ii) (also known as h 0 (n )) is that of a lowpass filter with its cutoff frequency at the 
Nyquist frequency rate ( H(n) = 0 ). It can also be shown that the frequency response of 
the wavelet coefficients ]\ ( n ) is that of a highpass filter with its maximum at 

H(n) = V2 and its zero magnitude at H{ 0) . This is the property that allows the input 
signal to be filtered into a “constant-Q” filter bank resulting in equal highpass and 
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lowpass portions preceding every downsample as described in Figure III-6. Figure III-S 
shows the frequency response of the PR filter bank used for the Haar Wavelet system. 

1. Theorem 5: 

If (pit) is a L l solution of the MRA equation (3.14), then the following equivalent 
conditions must be true: 




(A. 17) 



n 



For the Haar Wavelet the above equation leads to: 




h(0) + h(l) = H(0) = V2 . 



2. Theorem 6: 



For h(n)e L 1 , then 



h(2n) = h(2n + 1) if and only if H(n) = 0 . 



(A. 18) 



n 



n 



For the Haar Wavelet the above equation leads to: 
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3. 



Theorem 7: 



If <p(t) is a solution to the MRA equation (3.14) in L 2 nL 1 and is a 

solution of 



<S>(co) = -^H 

V2 



fco} 

9 

~ J 






fa\ 

yf’ 



(A. 19) 



and after iteration becomes 



&(*>) = n 

k = 1 





3>(0), 



(A. 20) 



then 



J (p(t)(p(t - k)dt = S(k) if and only if + =1. (3.21) 

1 

This theorem is the frequency domain equivalent to the time domain definition of 
the orthogonality of the scaling function [8]. Equation (3.42) shows that for a function 
that is narrow in time is spread out in frequency. 

For the Haar Wavelet the above equation leads to: 

<f>(a>) is the Fourier Transform of the scaling function (p(t) . Finding is an 

iterative process and equation (3.40) is often called the “Cascade Theorem”. By 
assuming an initial d>(0), the use of a computer to iterate the process can show that 
<t fco) converges, its inverse Fourier Transform is the scaling function and equation 
(3.42) is satisfied. As stated before, once the scaling function is determined, then by 
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equation (3.18) one can determine the shape of the Mother Wavelet. The shape of the 
Haar Scaling and corresponding Wavelet function is shown in Figure III-7. 



Theorem 8: 



For any h(n) E L 1 , 



^h(n)h(n-2k)=S(k) if and only if \H(o)f +\h(co + nf =2. (3.47) 



For the Haar Wavelet the above equation leads to: 

Starting with/f (a>) = -^=(l + e~ Jto ), we take the absolute value as multiplying by 
\ 2 

its complex conjugate giving, 



1 1 
■ + — -=e 






-JO) 



a/2 V2 U/2 42 



1 1 
■ + —=e 



JO) 



- + -e jw +-e' jm +- 
2 2 2 2 



Knowing that the complex exponential form of the cosine function is 



cos(x) = ±(e b! +e-“), 



we substitute this into the above expression giving 



|//(^)| 2 =l + cos(w). 



Now taking H ( co + tt) = -|= (l + e liu> * n) ) and multiplying by its complex 

V 2 



conjugate gives, 
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1 



1 



V2 H 



-j{(U+7t) 



1 



1 



^ 4i e 



j(C0+7t) 



1 J_ 

2 2 



= - + -e*V* + -<? 



9 




Again we substitute the cosine function in for its complex exponential form 

except now it is shifted by 180 degrees because of Euler’s Identity e” =-l. The 
resulting expression yields 



\H (a) + n)\ 2 = 1 - cos(ty) . 

Now we add both expressions together resulting in the final expression 
\H(co)\' +\h(co + 7zJ~ = 1 + cos(ty) + 1 -cos(ty) = 2. 
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APPENDIX B. FURTHER RESEARCH 



The denoising procedure used in Method 1 uses the Haar Wavelet system, which 
is the most basic wavelet set. Looking at the wavelet decomposition as an inner product 
between the basis set and the signal, we see that the Haar wavelet (Figure III-8) works 
well with functions that are piecewise continuous, but not smooth as in the case of most 
signals. Therefore it would be better to find a basis set where the Mother Wavelet has a 
shape similar to that of the signal, such as a wavelet from the Daubechies family [15] 
shown in Figure B-l. By denoising using a smoother, compact wavelet we still get the 
highpass and lowpass effect, but now the wavelet coefficients (output of the filter) might 
be larger due to better correlation between the signal and wavelet. Consider the Haar 
Wavelet in the frequency domain shown in Figure III-9. It has been shown that the Haar 
filter coefficients are only first-order polynomials. These basic filter coefficients imply 
that the cutoff frequencies of the highpass and lowpass channels overlap significantly. 
Therefore, the DWT denoising process might perform better with a filter pair whose 
frequency response is sharper at the cutoff frequency like Daubechies 10 (dblO) as shown 
in figure B-2. Indeed, upon one simple test using the same 5 kHz transient signal, the 
dblO wavelet system was able to reduce the noise threshold lower than the Haar wavelet, 
as shown in Figures B-3 through B-6. 
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Figure B-l Daubechies 10 (dblO) Wavelet 




Figure B-2 Frequency Response of dblO Wavelet System (f s =l 1.025 kHz) 
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Figure B-3 Wavelet Packet Framework ( AAD 3 and DAD 3 ) 
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Figure B-4 Denoising 5 kHz Signal Using Haar Wavelet System ( AAD 3 and DAD 3 ) 
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Figure B-5 Wavelet Packet Framework ( AAD 3 and DAD 2 ) 
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Figure B-6 Denoising 5 kHz Signal Using dblO Wavelet System ( AAD 3 and DAD 3 ) 
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Another modification to Method 1 that might be considered for further research is 
to define a Figure of Merit (FOM), which would allow for comparison of the various 
wavelet channels. To this end, we investigate the statistics of white noise. By definition, 
a zero mean white noise process x(?)is defined by its first moment (or mean) described 
as 



£{*(*)}= 0 (B.l) 

and second moment (or variance) described as 

E{x{t)x(t-T)}=G W 2 6{T). (B.2) 

In general, we wish to come up with an operation that can describe the correlation 
between N white noise processes, where N is the number of channels in the wavelet 
decomposition and the resulting correlation matrix R xx is of dimension N xN . In 

addition, this operation must provide visual clarity in distinguishing a transient signal 
from background noise. With this in mind, a figure of merit is developed as 



FOM =1- 



min eig(R xx ) 
max eig(R xx ) 



(B.3) 



The ratio of the minimum eigenvalue to the maximum eigenvalue will always be 
a number between zero and unity. If the minimum eigenvalue is equal to the maximum, 
indicating a set of completely uncorrelated signals with normalized power, then the ratio 
would equal 1 and the whole FOM would equal zero. On the other end, if there is some 
correlation between any two signals (i.e. a transient detected in at least one channel of 
the multiresolution filter bank), the ratio in the expression (B.3) will be less than unity 
and the entire expression will be greater than zero. This non-zero value will show up as a 
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“spike” when the FOM is plotted as a function of time. The square root simply smoothes 
out the FOM so that an indication of a transient signal is easier to see on a visual display. 

There is one final note about this method. Because of the downsampling 
operation after the filter bank, the length of data will decrease by a factor of two at each 
stage. In order to make every vector the same length, the signal is upsampled and run 
through an averaging filter. The effect of this process “stretches” the signal to the desired 
number of samples, and makes every other sample the average of its two neighbors. 
Performing this operation does not affect the statistics of the signal, nor does it wash out 
any transient information. The other option to perform the correlation operation is to 
downsample every detail to the size of the finest approximation. This too would be a 
valid process, but with downsampling comes the loss of information. The other reason to 
upsample and average is that the correlation function determined by a computer is an 
averaging operation. The upsample-and-average method provides more data points to 
average and therefore a better correlation representation. After calculating the details and 
approximations, the data is run through an upsample/averaging process to achieve the 
proper data length (indicated as primed details and approximations in Figure B-7). 

The re-sizing of the detail and approximation at the finest scale are then combined 
together as a set of column vectors as shown in Figure B-8 and in the equation 



t t 



t 



vv = 



i i 



(B.4) 



i 
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Figure B-7 Modified Method 1 Process 
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For this procedure, three stages of multiresolution are again used. Finally, the 
correlation function is computed as 



„ 1 T 

K w =~ w 



(B.5) 



and the figure of merit is determined as 



FOM(t) = 1- 



I min eig{R„Jtjj 
\ maxejg(/?^.(0) ' 



(B.6) 



A plot of the Figure of merit reveals values going from zero to unity as correlation 
between any two inputs increases. A “spike” on the FOM plot reveals a transient signal. 
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APPENDIX C. MATLAB PROGRAMS 
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A. 



CODE COMMON TO BOTH METHODS 



1. ARMA model 

function [b, a] =cristi_ar (y ) ; 

% 

% [b, a] =cristi_ar (y ) ; 

o 

% y = data (a COLUMN vector) 

ry=cris ti_corr (y , y ) ; 

% Levinson recursion 

a=l; b=l; sig2=ry(l); sig2p=sig2; 

p= 0 ; f ac t = 0 ; 

% while ( p< 10); 

while (p<10 ) & ( fact< 0 . 98 ) ; 

r=ry (2 :p+2) ; 

D=r' *flipud(a) ; Dprime=r ' *flipud(b) ; 
g=D/sig2p; gprime=Dprime/sig2 ; 

ap= [a; 0] -g* [0; flipud(b) ] ; b= [b; 0] -gprime* [0; flipud(a) ] 
a=ap; 

fact=l-g*gprime; 
sig2=fact*sig2 ; 
sig2p=fact *sig2p; 

P=P+1; 

end 

b=sqrt ( sig2p) ; 



2. Autocorrelation function 

function Rxy=crist i_corr (x,y) 

%r unction will calculate the cross-correlation function 
%It is faster than MAT LAB function XCORR 
% Rxy=cristi_corr (x # y) 

% 

m=max ( [ length (x) , length (y ) ] ) ; 
fx=f f t (x, 2 ^m) ; 
fy=f f t (y, 2*m) ; 

Sxy=fx. *conj (fy) ; 

Rxy=real ( if f t (Sxy) ) /m; 

%end % cristi_corr 
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3 



Benchmark Test 



%Name: John D. Stevens 

%Date: 22 Jul 1999 

%Filename : meth_0102 .m 

%description : this program uses Methods 1 and 2 to detect 

various 

%transient signals in a colored noise environment. 

a 

o 

%====================================================:============ 

clear 

% [s, fs,bps] =wavread( 'H: \matlab2\ thesis \data\wav_data\ record 7 ) ; 
%school 

% [S, fs,bps] =wavread( ' C: \My 

Documents \dad\m_f iles\ thesis\data\wav_data\record / ) ; %home 
% load ' H : \matlab2 \ thesis \ data \mat_data\ fan .mat' 
load ' H : \matlab2 \ thesis \ da ta\mat_data\ fan2 .mat ' 

%load ' C : \My Documents \cad\m_f iles\ thesis\data\wav_data\ f an . mat ' 
%load ' C : \My Documents \dad\m_f iles\ thesis\data\wav_data\ f an2 . mat ' 
fs=11025; % telephone quality 

%initialize the confidence vector: 

%narrowband signals: 
nbl= [50; .5; .5] ;nblh=[50; .1; .5] ; 
nb2=[200; .1; .1] ;nb2h=[200; .2; .2] ; 
nb3= [600; .1; .05] ;nb3h=[600; . 1; .08] ; 
nb4 =[900; .1; .03] ;nb4h=[900; . 1; .05] ; 
nb5= [2000; .1; .01] ;nb5h=[4000; .1; .1] ; 
nb6= [5000 ; .1; .01] ;nb6h=[5000; .1; .003] ; 

NB= [nbl , nb2 , nb3 , nb4 , nb5 , nb6 ] ; 

NBH= [nblh / nb2h / nb3h / nb4h / nb5h,nb6h] ; 

BBC= [ .1; .8] ; 

BBW= [ .1; .01] ; 
trials=length (NB) ; 
max_dtr=max(NB / ) ; 
max_dtr=f loor (fs*max_dtr (2) ) ; 

fig=l; 

%TRSOUND= [ ] ; 

%CSOUND= [ ] ; 

%WSOUND= [ ] ; 

%for c=l:trials; 
c = 6 ; 

%ftr=NB(l, c) ; 
ftr=NBH(l,c) ; 

%dtr=NB(2, c) ; 
dtr=NBH ( 2 # c) ; 

%dtr=B3C(l); %duration of colored noise transient 
%dtr=33W ( 1 ) ; %duration of white noise transient 

%atr=N3 ( 3 , c) ; 
atr=NBH ( 3 , c ) ; 
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%atr=33C ( 2 ) ; 
%atr=B3W(2) ; 



%create transient signal: 

%A) narrowband sinusoidal transient 

NS__T=f loor ( f s*dtr ) ; %number of samples in the transient 
tot_time=NS_T/f s ; 

t=0 : (tot_time/ (NS_T-1) ) :tot_time; 
y=cos (2*pi*ftr*t) ; 
ham_win=hamming (length (t) ) ; 
y=y . *ham_win 7 ; 

%3) broadband transient signal: 

%NS_?=f loor ( f s^dtr) ; ^number of samples in the transient 
%y=randn (1 ,NS_T) ; 

%ham_win=hamming (NS_T) ; 

%y=y . *ham_win 7 ; 

%C) chirp signal: 

%NS_T=f loor ( f s^dtr) ; % number of samples in the transient 
% tot_time=NS_T/ f s ; 

%t = 0 : ( tot_time/ (NS_T-1) ) : tot_time; 

%tp=[0 dtr/4 dtr/2 3^dtr/4 dtr] ; % time breakpoints 
%fp=[0 200 100 150 300] ; % instantaneous frequency 

breakpoints 

%p=polyf it ( tp, fp, 4) ; % fit 4th order polynomial over time 

%y=chirp ( t , p) ; 

%D) colored noise transient signal: 

%NS_T=f loor (f s*dtr) ; %number of samples in the transient 
%y=fan2 (20000 : 20000+NS_T-l } 7 ; 

%ham_wi n=hammi nc (NS_T) ; 

%y=y . *ham_win' ; 

%E) colored noise transient signal: ELVIS 

%load 7 C:\My 

Documents \dad\m_f iles\ thesis\data\wav_data\elvis_a .mat 7 
%load 7 C : \My 

Document s\dad\m_f iles\ thesis \data\wav_da ta\elvis_b. mat 7 
%load 7 H : \matlab2 \ thesis \data\mat_data\elvis_a .mat 7 
%load 7 H: \matlab2\thesis\data\mat_data\elvis_b .mat 7 
%NS_T=length(elvis_a) ; 

%y-elvis_a 7 ; 

%ham_win=hamming (NS_T) ; 

%y=y . *ham_win 7 ; 

%atr= . 3 ; 

yf = f an2 ; 

yf (10001: (10000+NSJT) )=yf (10001: (10000+NS_T) )+atr*y 7 ; 

% F) use the data fan (includes the tapping) as the input 

signal 

%yf =f an ; 

% met hod 1: pre-whitening 



152 



% [zO, zll_a, zlh_a, z21_a, z2h_a , z3 l_a , z3h_a j =method_01_f i 



(yf ) ; 



% [zO, zll_a, zlh_a, z21_a, z2h_a, z31_a, z3h_a, SNR_1] =method_01_f ilt 
v (yf ) ; 

% [zO, zll_a, zlh_a, z21_a, z 2 h_a , z31_a, z3h_a] =tree_filt (yf ) ; 

[ STGl , STG2 , STG3 , STG4 ] =wavepackz (yf); 

%SNR_l=my_SNRl (STG1 , STG2 , STG3 , STG4 ) ; 

%SNRl=my_SNR14plots ( STGl , STG2 , STG3 , STG4 ) ; 

%method 2: post-whitening (Cristi's algl) 

[ zll_b, zlh_b, z21_b, z2h_b, z31_b, z3h_b] =method_02_f il t (yf ) ; 



% [zll_b, zlh_b, z2 1_b , z2h_b, z3 l_b , z 3 h_b , S\ T R_2 ] =method_02_f iltv (y 
f ) ; 

%SNR2=my_SNR24plots (yf ) ; 

%method 3: SVD model (Cristi's alg2) 

% [zl_c, z2_c, z3_c, z4_c, z5_c, z6_c] =method_03_f lit (yf ) ; 



%y= [y, zeros (1, max_dtr- length (y) ) ] ; 

%TRSOUND= [ TRSOUND ; y ] ; 

%CSOUND= [CSOUND ; yf ' ] ; 

%WSOUND= [WSOUND; zO ' ] ; 

%plot results: 
s3_a=3 ; i3_a=3+s3_a; 
s3_b=4 ; i3_b=3+s3_b; 
s4_a=5 ; i4_a=7+s4_a; 
s4_b=6 ; i4_b=7+s4_b; 

f ig=plot_methO 102 ( STGl , STG2 ( 1 , :) / STG2(2 / :) , STG3(s3_a, :) , STG3(s 
3_b, : ) 

STG4 ( s4_a , : ) , STG4 (s4_b, : ) ,NS_T, fig) ; 

%f ig=plot_meth0102 { zO , zll_a , zlh_a , z21_a , z2h_a , z3 l_a , z 3 h_a , NS_T 

/ fig) ; 

f ig=plot_meth0102 (STGl, zll_b, zlh_b, z21_b, z2h_b, z31_b, z3h_b,NS_ 
T, fig) ; 

%f ig=plot_meth03 (zO, zl_c , z2_c, z3_c, z4_c , z5_c , z6_c , NS_T, fig) ; 
figure (fig) ; 

specgram (yf ,512, fs, 512, 400) ; 
f ig=f ig+1; 

%end 
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4. Plotting Routines 



%Name: John D. Stevens 

%Date: 26 Feb 2000 

%Fiiename: plot_meth0102 .m 

%description : this function will run a signal through the first 

algorithm 

%that was developed by Dr. Roberto Cristi and plot it 

a 

o 

% 

b=plot_meth0102 (zf0,zf0c,zfl,zflt,zf2 / zf2t,zfar, NS_T, f ig_num) ; 

% 

%================================================================ 

function 

b=plot_meth0102 (zfar, zf0,zf0t,zfl,zflt,zf2, zf2t / NS_T / f ig_num) ; 
TTL = ; 
slp=l ; 

%f ig_num=l ; 



% indexes due to downsampling : 

indl=10000 ; 

ind2=f loor ( indl/2 ) ; 

ind4=f loor ( indl/4 ) ; 

ind8=f loor ( indl/8 ) ; 

Ns_tl=NS_T; 

Ns_t2= floor (NS_T/ 2 ) ; 

Ns_t4=f loor (NS_T/4) ; 

Ns__t8 = f loor (NS_T/ 8 ) ; 

%create the spikes surrounding the signal 
sfar= zeros ( size ( zf ar) ) ; 
sfar (indl) =max (zfar) ; 

sfar ( indl + floor (Ns_tl*slp) ) =max (zfar) ; 

sf 0=zeros (size(zfO) ) ; 
s f 0 ( ind2 ) =max ( z f 0 ) ; 

sf 0 { ind2+ floor (Ns_t2*slp) ) =max ( zf 0 ) ; 

sf 0t=zeros (size (zf Ot) ) ; 
s f 0 1 { ind2 ) =max ( z f 0 1 ) ; 

sf Ot ( ind2 + f loor (Ns_t2*slp) )=max(zf0t) ; 

sf l=zeros (size(zfl) ) ; 
s f 1 ( ind4 ) =max ( z f 1 ) ; 

sf 1 ( ind4+f loor (Ns_t4* sip) ) =max ( zf 1 ) ; 

sf lt=zeros (size { zf it ) ) ; 
s f 1 1 ( ind4 ) =max ( z f 1 1 ) ; 

sf It ( ind4+f loor (Ns_t4* sip) ) =max { zf It ) ; 

sf 2= zeros (size(zf2) ) ; 
sf 2 { ind8 ) =max { zf 2 ) ; 
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sf 2 ( ind8+ floor (Ns_t8*slp) ) =max ( zf 2 ) ; 

sf2t=zeros (size(zf2t) ) ; 
sf2t ( ind8) =max (zf2t) ; 

sf 2 1 ( ind8+ floor (Ns_t8 *slp) ) =max(zf2t) ; 

^create the time vectors: 

Nar=length ( zfar) ; 

N0=length ( zf 0 ) / 

Nl=length ( zf 1) ; 

N2=length(zf2) ; 

fs=11025; 
f_time=Nar/f s; 

t_ar=0 : f_time/ (Nar-1) : f_time; 
t_0 = 0 : f_time/ (NO-1 ) : f_time; 

t 1 = 0 : f__time/ (Nl-1) : f__time; 

t_2 = 0 : f_time/ (N2-1) : f_time; 

%plot the figures 
figure ( f ig_num) 

plot ( t_ar, zfar, 'b' , t_ar, sfar, ' r: ' ) ; 
title ( [TTL, ' AR WHITENING FILTER ' ] ) ; 
ylabel ( ' MAGNITUDE # ) ; 

axis ( [0, f_time, min (zfar) # max (zfar) ] ) ; 
xlabel ( 'TIME (seconds) ' ) ; 

%ylabel ( ' MAGNITUDE ' ) ; 

% legend ( [ ' SNR = ' , num2str ( SNR( 1 ) ) , ' dB']); 
f igure ( f ig_num+ 1 ) 

subplot (2,1,1) , plot ( t_0 , zf 0 , 'b ' , t_0 , sf 0 , ' r : ' ) ; 

%title( [TTL, 'LOWPASS FILTER STAGE 
1 ' ] ) , axis ( [ 0 , f_time, min ( z f 0 ) , max ( zf 0 ) ] ) ; 
title ( [TTL, 'LOWPASS FILTER 
(APPROXIMATION) ' ] ) ,axis ( [0, f_time, 0,1]); 

%ylabel ( 'MAGNITUDE' ) ; 

% legend ([' SNR = ' , num2s tr ( SNR ( 2 ) ) , ' d3' ] ) ; 
subplot (2,1,2), plot ( t_0 , zf Ot , ' b' , t_0 , sf Ot , ' r : ' ) ; 

%title( 'HIGHPASS FILTER STAGE 
1 ' ) , axi s ( [ 0 , f _ t ime , min ( z f 0 1 ) , max ( z f 0 1 ) ] ) ; 

title ( [TTL, 'HIGHPASS FILTER ( DETAIL) ' ] ) , axis ( [ 0 , f_time, 0,1]) 
xlabel ( 'TIME (seconds) ' ) ; 

%ylabel ( 'MAGNITUDE' ) ; 

% legend ( [ ' SNR = ' , num2 s t r ( SNR ( 3 ) ) , ' d3']); 
figure (f ig_num+2) 

subplot (2,1,1) , plot ( t_l , zf 1 , ' b' , t_l , sf 1 , ' r : ' ) ; 

%title( [TTL, 'LOWPASS FILTER STAGE 
2 ' ] ) , axis ( [ 0 , f_time, min ( zf 1 ) , max ( zf 1) ] ) ; 
title ( [TTL, 'LOWPASS FILTER 
(APPROXIMATION) ' ] ) ,axis ( [0, f_time, 0, 1] ) ; 

%ylabel ( 'MAGNITUDE' ) ; 

% legend ( [ ' SNR = ' , num2str ( SNR ( 4 ) ) , ' d3']); 
subplot (2,1,2) , plot ( t_l , zf It , ' b' , t_l, sf 1 1 , ' r : ' ) ; 
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%ti tie ( ' HIGHPASS FILTER STAGE 
2 ' ) , axi s ( [ 0 , f _t ime , min ( z f 1 1 ) , max ( z f 1 1 ) ] ) ; 

title ( [TTL, 'HIGHPASS FILTER (DETAIL) ' ]) , axis ([ 0 , f_t ime , 0 , 1 ]) ; 
xlabel ( 'TIME (seconds) ' ) ; 

%ylabel ( ' MAGNITUDE ' ) ; 

%legend( [ ' SNR = ' , num2str ( SNR ( 5 ) ) , ' dB ' ] ) ; 
figure ( f ig_num+3 ) 

subplot (2,1,1) , plot ( t_2 , zf 2 , 'b ' , t_2 , sf 2 , ' r : ' ) ; 

%title( [TTL, 'LOWPASS FILTER STAGE 
3 ' ] ) ,axis ( [0, f_time,min ( zf 2 ) ,max(zf2) ] ) ; 
title ( [TTL, 'LOWPASS FILTER 
(APPROXIMATION) ' ] ) , axis ( [0, f_time, 0,1]) ; 

%ylabel ( ' MAGNITUDE ' ) ; 

% legend ( [ ' SNR = ' , num2str ( SNR ( 6 ) ) , ' dB']); 
subplot (2,1,2) , plot ( t_2 , zf 2t , ' b' , t_2 , sf 2 1 , ' r : ' ) ; 

%title( 'HIGHPASS FILTER STAGE 
3 ' ) , axis ( [ 0 , f_time , min ( zf 2t ) , max ( zf 2 t ) ] ) ; 

title ( [TTL, 'HIGHPASS FILTER (DETAIL) ']), axis ([ 0 , f_t ime , 0 , 1 ]) ; 
xlabel ( ' TIME ( seconds ) ' ) ; 

%ylabel ( 'MAGNITUDE' ) ; 

%Iegend ( [ ' SNR = ' , num2str ( SNR ( 7 ) ) , ' d3'l); 



b=f ig_num+4 ; 



%Name : John D. Stevens 

%Date : 22 Jul 1999 

%Filename: SNR_plots3 .m 

%description : plots the results from my_SNRl and my_SNR2 
%======================================================= 

clear 

%load up results: 

load 'H : \matlab2\ thesis \data\mat_data\SNR_ARb. mat ' 
load 'H: \matlab2 \ thesis \ data \mat__data\SNR_Mlb . mat ' 
load ' H : \matlab2 \ thesis \data\mat_data \SNR_M2b. mat ' 

%load 'C:\My 

Document s \ dad\m_f iles\thesis\data \wav_da t a \ snr_a r . ma t ' 
%load 'C:\My 

Document s \dad\m_f iles\ thesis \data\wav_data\snr_Ml .mat ' 
%load 'C:\My 

Documents \dad\m_f iles\ thesis\data\wav_data\snr_M2 .mat ' 

FTR= 0:10:5000; 

ATR=0 : . 002 : 1; 

thresh=3; %threshold in dB 

ax= [min ( FTR) , max ( FTR) ,min(ATR) , max ( ATR) ,-20,50]; 



figure ( 1 ) 

%surf (FTR, ATR, SNR_AR ) ;axis (ax); 
mesh ( FTR, ATR, SNR_AR, SNR_M1) ;axis(ax) ; 

C=colormap ; 

title ( 'NARROWBAND DETECTION REGION FOR WHITENING FILTER'); 
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xlabel( 'FREQUENCY OF TRANSIENT (Hz) ' ) ; 
ylabel ( 'AMPLITUDE OF TRANSIENT' ) ; 
z label ( ' SNR (d3) ' ) ; 

figure (2 ) 

%surf (FTR, ATR, SNR_M1) ; axis (ax) ; 
mesh (FTR, ATR, SNR_M1 , SNR_M1 ) ; axis (ax) ; 

title ( ' NARROW3AND DETECTION USING WAVELET DENOISING ' ) ; 
xlabel ( 'FREQUENCY OF TRANSIENT (Hz)'); 
ylabel ( 'AMPLITUDE OF TRANSIENT' ) ; 
z label ( ' SNR (dB) ' ) ; 

figure (3 ) 

% surf ( FTR , ATR , SNR_M2 ) ; axi s ( ax ) ; 
mesh ( FTR , ATR , SNR_M2 , SNR_M1 ) ; axis (ax); 

title ( 'NARROWBAND DETECTION USING MULTIRESOLUTION INNOVATION - 
THRESHOLD = 3 dB'); 

xlabel (' FREQUENCY OF TRANSIENT (Hz)'); 
ylabel ( 'AMPLITUDE OF TRANSIENT'); 
z label ( ' SNR (dB) ' ) ; 

figure (4) 

X=SNR_AR ; 

I=f ind (x<thresh) ; 
x(I) =thresh; 

contourf (FTR, ATR,x) ;colorbar; 
colormap( jet) ; 

%pcolor (FTR, ATR,x) ; 

title ( 'NARROWBAND DETECTION REGION FOR WHITENING FILTER - 
THRESHOLD = 3 d3 ' ) ; 
xlabel ( ' FREQUENCY (Hz ) ' ) ; 
ylabel ( 'AMPLITUDE' ) ; 

figure (5) 
x=SNR_Ml ; 

I=f ind (x<thresh) ; 
x(I) =thresh; 

contourf (FTR, ATR, x) ;colorbar; 
colormap ( jet) ; 

%pcolor (FTR, A.TR, X) ; 

title ( 'NARROWBAND DETECTION USING WAVELET DENOISING - THRESHOLD 
3 dB ' ) ; 

xlabel ( ' FREQUENCY ( Hz ) ' ) ; 
ylabel ( 'AMPLITUDE' ) ; 

figure ( 6 ) 
x=SNR_M2 ; 

I=f ind (xcthresh) ; 
x(I) =thresh; 

contourf (FTR, ATR, x) ; colorbar 
colormap (jet) ; 

%pcolor (FTR, ATR, X) ; 

title ( 'NARROW3AND DETECTION USING MULTIRESOLUTION INNOVATION - 
THRESHOLD = 3 c3 ' ) ; 
xlabel ( ' FREQUENCY (Hz ) ' ) ; 
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ylabel ( 7 AMPLITUDE ' ) ; 



%DETERMINE PERCENT IMPROVEMENT ON DETECTION: 

%A) FOR THE ENTIRE SPECTRUM: 
tot_pts_AR=prod (size ( SNR_AR) ) ; 

det__pts_AR=f ind ( SNR_AR>= thresh) ; %find all points >= 3db; 
det__pts_AR= length ( det__pts_AR) ; 

pc t_det_AR= ( det_.pt s_AR/ to t__pts_AR) *100; ^percent detection for 
whitening filter 

tot_pts_Ml=prod (size ( SNR__Ml ) ) ; 

det_pts_Ml = f ind(SNR_Ml> = thresh) ; %find all points >= 3db; 
det_pts_Ml=length (det_pts_Ml ) ; 

pct_det_Ml= (det_pts_Ml/tot_pts_Ml ) *100; %percent detection for 
whitening filter 

pct__incrMl=100* (pct_det_Ml-pct_det_AR) /pct_det_AR; %detection 
percent increase 

tot_pts__M2=prod(size ( SNR_M2 ) ) ; 

det_pts_M2=find(SNR_M2>=thresh) ; %find all points >= 3db; 
det_pts_M2=length (det_pts_M2 ) ; 

pct_det_M2= (det_pts_M2/tot_pts_M2 ) *100 ; %percent detection for 
whitening filter 

pct_incrM2=100* (pct_det_M2 -pct_det_AR) /pct_det_AR; %detection 
percent increase 

PCT_DET_TOT= [pct_det_AR; pct_det_Ml ;pct_det_M2] 

PCT_BETTER= [ 0 ; pct__incrMl ; pct_incrM2 ] 

%3) FOR LOW FREQUENCIES BELOW 2kHZ ONLY: 
f in=f ind (FTR==2000) ; %column where freq = 2000 Hz 

tot_2k_AR=prod (size (SNR_AR ( : # 1 : fin) ) ) ; 

det_2k_AR=f ind ( SNR_AR ( : , 1 : f in) >= thresh) ; %find all points >= 
3db; 

det_2k_AR=length (det_2k_AR) ; 

pct__2k__AR= (det_2k_AR/ tot_2k_AR) *100; %percent detection for 
whitening filter 

tot_2k_Ml=prod ( size (SNR_Ml ( : , 1 : fin) ) ) ; 

det_2k_Ml=f ind ( SNR_M1 ( : , 1 : f in) >=thresh) ; %find all points >= 
3db; 

det_2k_Ml=length (det_2k_Ml ) ; 

pct_2k__Ml = (det_2k__Ml/tot_2k_Ml ) *100 ; %percenz detection for 
whitening filter 

pct_2k_incrMl=100* (pct_2k_Ml-pct_2k__AR) /pct_2k_AR; %detection 
percent increase 

tot_2k_M2=prod ( size ( SNR_M2 ( : , 1 : f in) ) ) ; 

det_2k_M2=f ind ( SNR_M2 ( : , 1 : f in) >=thresh) ; %find all points >= 
3db ; 

det_2k_M2= length (det_2k_M2 ) ; 

pct_2k__M2 = (det_2k_M2/ tot_2k_M2 ) *100 ; %percent detection for 
•whitening filter 
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pct_2k_incrM2=100* (pc t_2k_M2 -pct_2k_AR) /pct_2k_AR; %detection 
percent increase 

PCT_DET_2 k= [ pc t_2 k_AR ; pc t_2 k_Ml ; pc t„2 k_M2 ] 

PCT_BETTER_2k= [ 0 ; pct_2k_incrMl ;pct_2k_incrM2 ] 



%C) FOR LOW FREQUENCIES BELOW 300 Hz ONLY: 
f in=f ind ( FTR==300 ) ; %column where freq = 300 Hz 

tot_3c_AR=prod ( size (SNR_AR( : , 1 : fin) ) ) ; 

det_3c_AR=f ind ( SNR_AR ( : , 1 : f in) >= thresh) ; %find all points >= 
3db; 

det_3c_AR= length (det_3c_AR) ; 

pct_3c_AR= (det_3c_AR/ tot_3c_AR) *100 ; %percent detection for 
whitening filter 

tot_3c_Ml=prod(size(SNR_Ml ( : , 1: fin) ) ) ; 

det_3c_Ml=find(SNR_Ml ( : , 1 : fin) >=thresh) ; %find all points >= 
3db ; 

det_3c_Ml=length(det_3c_Ml) ; 

pct_3c_Ml= (det_3c_Ml/tot_3c_Ml) *100 ; ^percent detection for 
whitening filter 

pct_3c_incrMl=100* (pct_3c_Ml-pct_3c_AR) /pct_3c_AR; %detection 
percent increase 

tot_3c_M2=prod ( size ( SNR_M2 ( : , 1 : fin) ) ) ; 

det_3c_M2=f ind ( SNR_M2 ( : , 1 : f in) >=thresh) ; %find ail points >= 
3db ; 

det_3c_M2=length(det_3c_M2 ) ; 

pct_3c_M2 - (det_3c_M2/ tot_3c_M2 ) *100 ; ^percent detection for 
whitening filter 

pct_3c_incrM2=100 * (pct_3c_M2-pct_3c_AR) /pct_3c_AR; %detection 
percent increase 

PCT_DET_3c= [pct_3c_AR;pct_3c_Ml;pct_3c_M2] 

PCT_BETTER_3 c= [ 0 ; pc t_3 c_incrMl ; pc t_3 c_incrM2 ] 
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3 . 



Display Filter 



%Name: John D. S t evens 

%Date : 21 Jan 2000 

%Fiiename : smoother .m 

%description : this program will take a white noise process and 

output 

%a normalized smoothed version for display. This process also 
%gives an approximation to the noise variance and the signal 
variance, 

%which is useful to compute the signal to noise ratio 

% 

% zf =smoother (y) 

% 

%=============================================================== 

function zf =smoother (y ) 

if size (y, 2 ) >size (y , 1 ) , y=y' ; end % make y into column vectors 

Z=y. "2; 

B=1 ; 

A= [ 1 , - . 99 ] ; 

zf= filter ( 1 , [1, — .99] ,z) ; 
z f = z f ' ; %make output ir. row vectors 
%normalize the results: 
mzf =max (zf , [ ] ,2) ; 

zf=zf . /mzf ( : , ones (length(zf ) , 1) ) ; 

% [ z f , f_f in] =f ilter (3, A, z , y_init ) ; 



6. Broadband Transient SNR Test 

%Name : John D. Stevens 

%Date : 22 Jul 1999 

%Filename : SNR_elvis 

%description : 

%================================================================ 

clear 

load 'H: \matlab2\ thesis\data\mat_data\fan2 .mat' 
load ' H : \mat lab2 \ thesis \data\mat_data\ elvis_a . mat ' 

%load 'C : \My Documents \dad\m_f iles\ thesis\data\wav_data\fan2 .mat ' 
%load 'C:\My 

Documents \dad\m_f iles\thesis\data\wav_data\elvis_a . mat ' 

fs=11025; % telephone quality 

SNR_AR= [ ] ; 

SNR_M1= [ ] ; 

SNR_M2= [ ] ; 

dtr= . 1 
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%NS_T=length ( elvis_a) ; 

NS_T=f loor ( f s*dtr) ; % number of samples in the transien 

ham_win=hamming (NS_T) ; 

% ELVIS 

%y=elvis_a ' ; 

%y=y/max(y); %normalize the eivis signal 
%y=v . ~ham_win ' ; 

% WHITE NOISE: 
y=randn ( 1 , NS_T ) ; 
y=y . *ham_win' ; 

for atr=0 : . 005 : 1 ; %amplitude 
atr 

%add transient to benchmark signal 
yf =f an2 ; 

yf (10 001 : ( 10 00 0+NS_T) ) =yf (10 001 : (10000+NS_T) ) +atr*y ' 

%mechod 1: pre-whitening 
[ STG1 , STG2 , STG3 , STG4 ] =wavepack (yf ) ; 
snr_l=my_SNRl ( STG1 , STG2 , STG3 , STG4 ) ; 

%method 2: post-whitening 
snr_2 =my_SNR2 (yf ) ; 

SNR_AR= [ SNR_AR ; snr_l ( 1 ) ] ; 

SNR_M1= [ SNR_M1 ; snr_l ( 2 ) ] ; 

SNR_M2 = [ SNR_M2 ; snr_2 ] ; 

end 

ATR= (0 : . 005 : 1) ' ; 

res = [ SNR_AR , SNR_Ml , SNR_M2 ] ; 

figure ( 16 ) 

plot (ATR, res , ATR, 3*ones (size (ATR) ) , 'k: ' ) ; 
axis([0 1 -20 max ( SNR_M1 ) ] ) ; 
xlabel ( 'AMPLITUDE OF ELVIS'); 
ylabel ( ' SNR (dB)'); 

legend ([ 'Whitening Filter '],[' Method 1'], ['Method 2']); 
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B. METHOD 1 



1. Method 1 Filter Process 

%Name: John D. Stevens 

%Date: 1 March 2000 

%Filename: wavepack.m 

%description : this program will first whiten a colored noise 
signal 

%y and then decompose it into a Haar wavelet decomposition packet 
% framework the channels will be white noise processes 
% 

%================================================================ 

function [ STAGE 1 , STAGE 2 , STAGE 3 , STAGE 4 ] =wavepack (y ) 

if size (y, 2 ) >size (y, 1) , y=y ' ; end %make the signal a column vector 

%perform ARMA process to get a whitened signal 

[b, a] =cristi_ar (y) ; 

e0=filter (a,b,y) ' ; 

%f liter the signal through the Haar wavelet transform packet: 

M= [1/sqrt ( 2 ) , 1/sqrt (2 ) ; 1/1/sqrt ( 2 ) , -1/sqrt (2 ) ] ; 
temp=M*reshape (eO , 2 , length (e0) /2 ) ; 

Al=temp ( 1, : ) ; 

Dl=temp (2 , : ) ; 

temp=M*reshape ( A1 , 2 , length ( A1 ) / 2 ) ; 

AA2=temp (1 , : ) ; 

DA2=temp (2 , : ) ; 

temp=M* reshape ( AA2 , 2 , length (AA2 ) /2 ) ; 

AAA3=temp (1, :) ; 

DAA3=temp(2 , ; 

temp=M*reshape (DA2 , 2 , length (DA2 ) / 2 ) ; 

ADA3=temp (1,:); 

DDA3=temp(2, : ) ; 

temp=M* reshape (D1 , 2 , length (Dl) / 2 ) ; 

AD2=temp(l / : ) ; 

DD2 = temp (2 , : ) ; 

temp=M* reshape ( AD2 , 2 , length ( AD2 ) /2 ) ; 

AAD3=temp ( 1 , : ) ; 

DAD3=temp (2 , : ) ; 

temp=M*reshape ( DD2 , 2 , length {DD2 ) / 2 ) ; 

ADD3=temp ( 1 , : ) ; 

DDD3=temp ( 2 , : ) ; 

STAGEl=eO ; 

STAGE2=[A1;D1] ; 

STAGE3= [AA2 ; DA2 ; AD2 ; DD2 ] ; 

STAGE4= [ AAA3 ; DAA3 ; ADA 3 ; DDA3 ; AAD3 ; DAD 3 ; ADD3 ; DDD3 ] ; 
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%Name: John D. Stevens 

%Date: 1 March 2000 

%Filename : wavepackz .m 

%description : this program will first whiten a colored noise 
signal 

%y and then decompose it into a Haar wavelet decomposition packet 
^framework. The result will include the display filtering 
process, 

%smoother 

% 

%================================================================ 

function [ STAGE 1 , STAGE2 , STAGE3 , STAGE4] =wavepack (y) 

if size (y , 2 ) >size (y, 1) , y=y ' ; end %make the signal a column vector 

%perform ARMA process to get a whitened signal 

[b, a] =cristi_ar (y ) ; 

e0=f ilter (a, b,y) ' ; 

%f liter the signal through the Haar wavelet transform packet: 

M= [1/sqrt ( 2 ) , 1/sqrt (2) ; 1/1/sqrt ( 2 ) , -1/sqrt (2 ) ] ; 
temp=M*reshape(eO, 2 , length (eO) / 2 ) ; 

Al=temp ( 1, : ) ; 

Dl = temp ( 2 , : ) ; 

temp=M*reshape ( Al , 2 , length (Al ) 12 ) ; 

AA2=temp ( 1 , : ) ; 

DA2=temp (2,:); 

temp=M*reshape (AA2 # 2 , length (AA2) /2) ; 

AAA3 = t emp ( 1 , : ) ; 

DAA3=temp (2 , : ) ; 

temp=M*reshape ( DA2 # 2 , length ( DA2 ) / 2 ) ; 

ADA3=temp ( 1 # : ) ; 

DDA3=temp (2 , : ) ; 

temp=M*reshape (Dl , 2 , length (D1 ) /2 ) ; 

AD2=temp ( 1 , : ) ; 

DD2=temp (2 , : ) ; 

temp=M* reshape ( AD2 , 2 # length (AD2 ) /2 ) ; 

AAD3=temp ( 1 , : ) ; 

DAD3=temp ( 2 , : ) ; 

temp=M*reshape (DD2 , 2 , length (DD2 ) / 2 ) ; 

ADD3=temp ( 1 , : ) ; 

DDD3 = temp(2 / : ) ; 

STAG21= smoother (eO) ; 

STAGE2=smoother ( [Al;Dl] ) ; 

STAGE3=smoother ( [ AA2 ; DA2 ; AD2 ; DD2 ] ) ; 

STAGE4=smoother ( [ AAA3 ; DAA3 ; ADA3 ; DDA3 ; AAD3 ; DAD 3 ; ADD3 ; DDD3 ] ) ; 
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2 , 



Narrowband Transient SNR Test 



%Name: John D. Stevens 

%Date : 11 Jan 2000 

%Filename : my_SNRl . m 

%descript ion : find the signal to noise ratio of the all the 

channels 

%in the Haar Wavelet Packet 

a 

o 

% SNR=my_SNRl ( STAGS1 , STAGE2 , STAGE3 , STAGE 4 } ; 

% 

%============================================================== 

function SNR=my_SNRl ( STAGE 1 , STAGE 2 , STAGE3 , STAGE4 ) ; 

%the assumption is that the transient signal is .1 seconds long 

NS_T=f loor ( .1*11025) ; 

slp=l; 

Ns_tl=NS_T; 

Ns_t2=f loor (NS_T/2) ; 

Ns_t4=floor (NS_T/4) ; 

Ns__t8 = f loor (NS_T/ 8 ) ; 

%zfar is a column vector while all the other inputs are row 
vectors 

% indexes due to downsampling: 
indl=10000 ; 

indl_f in=indl+f loor (Ns_tl*slp) ; 

ind2=f loor ( indl/2 ) ; 

ind2_f in=ind2+f loor (Ns_t2*slp) ; 

ind4=f loor ( indl/4 ) ; 

ind4_f in=ind4+f loor (Ns_t4*slp) ; 

ind8=floor (indl/8) ; 

ind8_f in=ind8+ floor (Ns_t8*slp) ; 

S= STAGE 1 ; 

A1=STAGE2 { 1 , : ) ; 

D1=STAGE2 (2, : ) ; 

AA2=STAGE3 ( 1 # : ) ; 

DA2=STAGE3 (2, : ) ; 

AD2=STAGE3 (3 , : ) ; 

DD2=STAGE3 (4, : ) ; 

AAA3=STAGE4 ( 1 , : ) ; 

DAA3=STAGE4 ( 2 , : ) ; 

ADA3=STAGE4 (3 , : ) ; 

DDA3 =STAGE4 (4 , : ) ; 

AAD3 =STAGE4 ( 5 , : ) ; 

DAD3=STAGE4 ( 6 , : ) ; 

ADD3=STAGE4 { 7 , : ) ; 

DDD3=STAGE4 { 8 , : ) ; 
k=l .5; 

nS=var { [ S ( 1 : ( indl -1 ) ) , S ( floor (k* indl_f in) : length ( S) ) ] ) ; 
nAl=var ( [ A1 { 1 : ( ind2 -1 ) ) , A1 ( floor (k* ind2_f in) : length (Al ) ) ] ) ; 
nDl=var ( [D1 ( 1 : ( ind2 -1 ) ) , D1 ( floor (k* ind2_f in) : length (D1 ) ) ] ) ; 
nAA2=var { [ AA2 ( 1 : ( ind4 -1 ) ) , AA2 { floor (k*ind4__f in) : length { AA2 ) ) ] ) ; 
nDA2=var ( [DA2 ( 1 : ( ind4-l ) ) , DA2 ( floor (k*ind4_f in) : length (DA2 ) ) ] ) ; 
nAD2=var ( [AD 2 ( 1 : ( ind4-l ) ) , AD2 (floor (k* ind4_f in) : length { AD2 ) ) ] ) ; 
nDD2=var ( [ DD2 ( 1 : ( ind4-l ) ) , DD2 (floor ( k* ind4_f in) : length (DD2 ) ) ] ) ; 
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nAAA3=var ( [AAA 3 ( 1 : ( ind8- 

1) ) , AAA3 ( floor ( k* ind8_f in) : length(AAA3 ) ) ] ) ; 
nDAA3 =var ( [DAA3 (1: (ind8- 

1) ) , DAA3 ( floor (k*ind8_f in) : lengch(DAA3 ))]); 
nADA3=var ( [ ADA3 (1: (ind8- 

1) ) , ADA3 (floor (k*ind8_f in) : length (ADA3 ) ) ] ) ; 
nDDA3 =var ( [ DDA3 ( 1 : ( ind8 - 

1) ) , DDA3 ( floor (k*ind8_f in) : length (DDA3 ))]); 
nAAD3=var ( [ AAD3 (1: (ind8- 

1) ) , AAD3 ( floor (k*ind8_f in) : length (AAD3 ) ) ] ) ; 
nDAD3 = var ( [ DAD 3 ( 1 : ( ind8 - 

1) ) , DAD3 ( floor (k*ind8_f in) : length (DAD3 ) ) ] ) ; 
n ADD 3= var ( [ ADD3 ( 1 : ( ind8- 

1) ) , ADD3 ( floor (k*ind8_f in) : length (ADD3 ) ) ] ) ; 
nDDD3 =var ( [ DDD3 ( 1 : ( ina8 - 

1) ) , DDD3 ( floor (k*ind8_f in) : length (DDD3 ) ) ] ) ; 

%FIND THE SIGNAL POWER: 
sS=var (S (indl : indl_f in) ) ; 
sAl=var (A1 ( ind2 : ind2_f in) ) ; 
sDl=var (D1 ( ind2 : ind2_f in) ) ; 
sAA2=var (AA2 ( ind4 : md4_f in) ) ; 
sDA2 = var ( DA2 ( ind4 : ind4_f in)); 
sAD2=var ( AD2 ( ind4 : ind4_f in) ) ; 
sDD2=var ( DD2 ( ind4 : ind4_f in) ) ; 
sAAA3=var (AAA3 { ind8 : ind8_f in) ) ; 

SDAA3 =var ( DAA3 ( ind8 : ina8_f in ) ) ; 
sADA3 =var ( ADA3 { ind8 : ind8_f in ) ) ; 

SDDA3 =var ( DDA3 ( ind8 : ind8_f in)); 
sAAD3=var (AAD3 ( ind8 : ind8_f in) ) ; 
sDAD3 =var ( DAD3 ( ind8 : ind8_f in) ) ; 
sADD3=var (ADD3 ( ind8 : ind8_f in) ) ; 
sDDD3 =var ( DDD3 ( ind8 : ind8_f in) ) ; 

%get the raw ratio: 
snr_S=sS/nS; 
snr_Al=sAl/nAl ; 
snr_Dl = sDl/nDl ; 
snr_AA2 = s AA2 / nAA2 ; 
snr_DA2 = sDA2 / nDA2 ; 
snr_AD2=sAD2/nAD2 ; 
snr_DD2 = sDD2/nDD2 ; 
snr_AAA3 = s AAA3 / nAAA3 ; 
snr_DAA3 = sDAA3/nDAA3 ; 
snr_ADA3 = s ADA3 / nADA3 ; 
snr_DDA3 =sDDA3 /nDDA3 ; 
snr_AAD3 = sAAD3 /nAAD3 ; 
snr_DAD3=sDAD3/nDAD3; 
snr_ADD3 =sADD3 /nADD3 ; 
snr_DDD3 =sDDD3 /nDDD3 ; 

C= [ snr _S ; snr_Al ; snr_Dl ; snr_AA2 ; snr_DA2 ; snr_AD2 ; snr_DD2 ; . . . 
snr_AAA3 ; snr_DAA3 ; snr_ADA3 ; snr_DDA3 ; snr_AAD3 ; snr_DAD3 ; snr_ADD3 ; sn 
r_DDD3 ] ; 

x=c;x=x-ones (size (x) ) ; 

I=f ind (x< . 01) ; 

x(I) = . 01 ; SNR=10 * loglO (x) ; 

SNR= [ SNR ( 1 ) ; max ( SNR ( 2 : 15) ) ] ; 
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c. 



METHOD 2 



1. Method 2 Filter Process 

%Name: John D. Stevens 

%Date: 11 Jan 2000 

%Fiiename : method_02_f ilt ,m. 

%description: this function will run a signal through the first 

algorithm 

%that was developed by Dr. Roberto Cristi 
% 

% [zl,zlt,z2,z2t,z3,z3t] =me thod_02_f ilt (y ) 

% 

%=============================================================== 

function [zl,zlt,z2,z2t,z3,z3t] =method_02_f ilt (y ) 

if size (y , 1 ) >size (y , 2 ) , y=y # ; end % y row vector 
k=l/sqrt (2) ; 

%k= . 5 ; 

[b, a] =cristi_ar (y ' ) ; 

[bl , al , bit , alt] =model2k (b, a, k) ; 

[b2 , a2 , b2t , a2t ] =model2k (bl , al , k) ; 

[b3 , a3 , b3t , a3t ] =model2k (b2 , a2 , k) ; 

M=k* [1,1; 1,-1] ; 

yl=M*reshape (y, 2 , length(y) / 2) ; 
y2=M*reshape (yl ( 1 , : ) , 2 , length (yl) /2 ) ; 
y3=M*reshape (y2 ( 1 # : ) , 2 , length (y2 )/2); 

%e0=f liter (a, b, y) ; 

el=f ilter (al # bl,yl(l,:)); 

elt=f ilter (alt , bit , yl (2 , : ) ) ; %row vector 

e2=f ilter (a2 # b2 / y2 (1,:)); 

e2t=f ilter (a2t , b2 t ,y2 ( 2 , : ) ) ; 

e3=f ilter ( a3 , b3 , y3 ( 1 , : ) ) ; 

e3t=f ilter (a3t , b3t , y3 ( 2 , : ) ) ; 

%z0= smoother (eO ) ; 
zl=smoother (el) ; 
zl=zl/max (zl) ; 
zlt=smoother (elt ) ; 
zlt=zlt/max ( zlt ) ; 
z2=smoother (e2 ) ; 
z2=z2/max ( z2) ; 
z2t=smoother (e2t ) ; 
z2t = z2t/ max ( z 2 1 ) ; 
z3=smoother (e3 ) ; 
z3=z3/max ( z3 ) ; 
z3t=smoother (e3t ) ; 
z3t=z3t/max(z3t); 
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2. Multiresolution Innovation Filter Bank Routine 

%Name: John D. Stevens 

%Date : 23 Feb 2000 

%Filename : model 2k . m 

%description : This program applies the multiresolution 

innovation process 
% 



function [b2 , a2 , b2t , a2t ] = model 2k (b, a, k) 

% [b2 , a2 , b2 t , a2t ] =moael2 (b, a) 

% "t" stands for ‘'tilde", which is the high pass filter 
[Ce, Co, Fz] =polyd (b, a) ; 

if length (Ce) < length (Co) , num=k*Co+k* [Ce, zeros ( 1 , length (Co) - 
length (Ce) ) ] ; end 

if length (Ce) >length (Co) , num=k*Ce+k* [Co , zeros ( 1 , length (Ce) - 
length (Co) ) ] ; end 

if length (Ce) ==length (Co) , num=k*Ce+k*Co ; end 
den=Fz ; 

if length (num) >length (Fz ) , den= [Fz , zeros ( 1 , length (num) - 
length (Fz) ) ] ; end 

[FI , G1 , Cl , Dl] =t f 2ss (num, den) ; 

if length(Co) +l<length (Ce) , num=k*Ce+k* [0,Co, zeros (1, length(Ce) - 
length (Co) -1) ] ; end 
if length (Co) +1> length (Ce) , 

num=k* [ 0, Co] +k* [Ce, zeros (1, length (Co) +1 -length (Ce) ) ] ; end 
if length (Co) +l==length (Ce) , num=k*Ce+k* [0,Co] ; end 
den=Fz ; 

if length (num) >length (Fz) , den= [den, zeros ( 1 , length (num) - 
length (den) ) ] ; end 

if length(num)<length(Fz) , num= [num, zeros (1, length(den) - 
length ( num) )] ; end 

[F2 , G2 , C2 , D2 ] =tf 2ss (num, den) ; 

% overall state space 
if length (F2) ==length(Fl) +1, 

Fl= [ [FI ; zeros ( 1 , length ( Fl ) -1 ) , 1 ] , zeros ( length ( FI ) +1 , 1 ) ] ; 

Gl= [G1 ; 0] ; 

C1=[C1, 0] ; 

end 



A=F1 ' ; 

B=[C2 # ,C1 # ] ; 

C=G2 ' ; 

D= [ D2 , Dl ] ; 

Q=B*B ' ; 

R=D*D ' ; 

S=B*D ' ; 

[K, P] = mydlqe (A, eye (size (A) ) , C, Q, R, S) ; 
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[b2,a2]=ss2tf(A,K,C,D ; 
sigma=sqrt (C*P*C ' +R) ; 
b2=b2/sigma ; 

% "tilde" - the high pass filtering 

if length (Ce) < length (Co) , num=k*Co-k* [Ce, zeros ( 1 , length (Co) - 
length (Ce) ) ] ; end 

if length(Ce) >length(Co) , num=k*Ce-k* [Co, zeros (1, length(Ce) - 
length (Co) ) ] ; end 

if length (Ce) ==length (Co) , num=k*Ce-k*Co ; end 
den=Fz ; 

if length (num) >length ( Fz ) # den= [Fz , zeros ( 1 , length (num) - 
length (Fz) ) ] ; end 

[FI , G1 , Cl , D1 ] =tf2ss (num, den) ; 

if length(Co) +l<length (Ce) , num=k*Ce-k* [0,Co, zeros (1, length(Ce) - 
length (Co) -1) ] ; end 

if length (Co) + 1> length (Ce) , num=k* [0,Co] - 
k* [Ce, zeros (1, length (Co) +1 - length (Ce) ) ] ; end 
if length(Co) +l==length (Ce) , num=k*Ce-k* [ 0 , Co ] ; end 
den=Fz ; 

if length (num) >length (Fz ) , den= [den, zeros ( 1 , length (num) - 
length(den) ) ] ; end 

if length (num) <length (Fz ) , num= [num, zeros ( 1 , length (den) - 
length (num) ) ] ; end 

[F2 , G2 , C2 , D2] =tf 2ss (num, den) ; 

% overall state space 
if length (F2 ) = = length (Fl) +1 , 

Fl= [ [FI; zeros ( 1 , length (Fl) -1 ) , 1] , zeros ( length (Fl) +1,1)]; 

Gl= [G1 ; 0 ] ; 

Cl= [Cl , 0 ] ; 

end 



A=F1 ' ; 

B= [C2 ' , Cl ' ] ; 

C=G2 ' ; 

D= [D2,D1] ; 

Q=B*B ' ; 

R=D*D ' ; 

S=B*D ' ; 

[K, P] = mydlqe (A, eye (size(A) ) , C , Q, R, S) ; 

[b2t,a2t]=ss2tf(A,K,C,l) ; 
sigma=sqrt (C*P*C 7 +R) ; 
b2 1 = b2 1 / sigma ; 
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3 



Narrowband Transient SNR Test 



%Name : John D. Stevens 

%Date : II Jan 2000 

% F i 1 ename : my_SNR2 . m 

%description : this function will determine the SNR of various 

%transients in white noise 
% 

% SNR=my_SNR2 (y) 

% 

%============================================================= 

function SNR=my_SNR2 (y) 

if size (y , 1 ) >size (y , 2 ) , y=y # ; end % y row vector 
kk=l/sqrt ( 2 ) ; 

%kk= . 5 ; 

[b, a] =cristi__ar (y ' ) ; 

[bl , al , bit , alt ] =model2k (b, a, kk) ; 

[b2 , a2 , b2 1 , a2 1 ] =model2k (bl , al # kk) ; 

[b3 , a3 , b3t , a3 t ] =model2k(b2 , a2 , kk) ; 

M=kk*[l,l;l,-1]; 
yl=M*reshape (y , 2 , length (y) 12 ) ; 
y2=M*reshape (yl ( 1 , : ) ,2 , length (yl ) 12 ) ; 
y3=M*reshape (y2 ( 1 # : ) ,2 , length (y2 )/2); 

%e0=f liter (a, b, y) ; 
el=filter (al,bl,yl(l, :)); 

elt = f ilter (al t , bit # yl ( 2 , : ) ) ; %row vector 

e2=f ilter (a2 , b2 , y2 ( 1 , : ) ) ; 

e2t=f ilter (a2t,b2t,y2 (2,:)); 

e3=f ilter (a3 , b3 ,y3 (1,:)); 

e3t = f ilter (a3t , b3t , y3 (2 # : ) ) ; 



%f ind SNR: 

NS_T=f loor ( . 1*11025 ) ; 
slp=l ; 

Ns_t 1=NS_T; 

Ns__t2 = f loor (NS_T/2 ) ; 

Ns_t4=f loor (NS_T/4 ) ; 

Ns_t 8=f loor (NS_T/8 ) ; 

%zfar is a column vector while all the other inputs are row 
vectors 

%indexes due to downsampling: 
indl=10000 ; 

indl_fin=indl + f loor (Ns_tl * sip) ; 
ind2=f loor ( indl/2 ) ; 
ind2_f in=ind2+f loor (Ns_t2 *slp) ; 
ind4=floor (inal/4) ; 
ind4_fm=ind4 + floor (Ns_t4*slp) ; 
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ind8=f loor ( indl/8 ) ; 

ind8_f in=ind8 + floor (Ns_t8 * sip ) ; 

zf 0=el ; 
zf Ot=elt ; 
zf l=e2 ; 
zf lt=e2t ; 
zf 2=e3 ; 
zf2t=e3 t ; 

k=l . 5 ; 

%find the variance of the noise 

mnO=var( [zf 0 (1: (ind2-l) ) , zfO ( floor (k*ind2_f in) : length (zfO) ) ] ) 
mnOt=var( [zfOt(l: (ind2- 

1) ) , zf Ot ( f loor (k*ind2_f in) : length (zf Ot ) ) ] ) ; 

mnl=var ( [zf 1 ( 1 : ( ind4-l ) ) , zf 1 (floor (k*ind4_f in) : length ( zf 1 ) ) ] ) 
mnlt=var ( [zflt (1 : (ind4- 

1 ) ) , zf It ( f loor (k*ind4_f in) : length ( zf It ) ) ] ) ; 

mn2=var ( [ zf 2 (1 : ( ind8-l) ) , zf2 (floor (k*ind8_f in) : length ( zf 2 ) ) ] ) 
mn2t=var ( [ zf 2 t ( 1 : ( ind8- 

1 ) ) , zf 2 t ( floor (k* ind8_f in) : length ( zf 2 t ) ) ] ) ; 

%find the variance of the spike in the region: 

mpO=var (zf 0 (ind2 : ind2_fin) ) ; 

mpOt=var ( z f Ot ( ind2 : ind2__f in) ) ; 

mpl=var (zfl (ind4 : ind4_fin) ) ; 

mplt=var ( zf It ( ind4 : ind4_f in) ) ; 

mp2=var (zf2 (ind8 : ind8_fin) ) ; 

mp2t=var ( zf 2t (ind8 : ind8_fin) ) ; 

%get the raw ratio: 
snr_0=mp0/mn0 ; 
snr_0t=mp0t /mnOt ; 
snr_l=mpl /mnl ; 
snr_lt=mplt/mnlt; 
snr_2=mp2 /mn2 ; 
snr_2t=mp2t /mn2t ; 



%SNR= [ snr_ar, snr_0 , snr_0t , snr_l , snr_lt , snr_2 , snr_2 t ] ; 
SNR= [ snr_0 , snr_0t , snr_l , snr_l t # snr_2 , snr_2 1 ] ; 

% check for ill conditions 
x=SNR; 

x=x-ones ( size (x) ) ; 

I=f ind (x< . 01 ) ; 
x(I) =.01; 

SNR=10 * loglO (x) ; 

SNR=max ( SNR) ; 



170 



LIST OF REFERENCES 



[1] K. L. Fracks, Improving Transient Signal Synthesis Through Noise Modeling 
and Noise Removal, Master’s Thesis, Naval Postgraduate School, Monterey, 
California, March 1994. 

[2] R. J. Urick, Ambient Noise in the Sea, Peninsula Publishing, Los Altos, 
California, 1986. 

[3] R. J. Urick, Principles of Underwater Sound , 3 rd Edition, McGraw Hill, Inc.: 
New York, 1983. 

[4] R. Leon-Garcia, Probability and Random Processes, Addison Wesley 
Publishing Company: Menlo Park, CA, 1989. 

[5] C. W. Therrien, Discrete Random Signals And Statistical Signal Processing, 
Prentice Hall: Englewood Cliffs, NJ, 1992. 

[6] L. Cohen, “Time - Frequency Distributions - A review,” Proceedings of the 
IEEE, Vol. 77, No. 7, pp. 941 - 981, July 1989. 

[7] O. Rioul and M. Vetterli, “Wavelets and Signal Processing,” IEEE Signal 
Processing Magazine, pp. 14 - 38, October 1991. 

[8] J. S. Lim and A. V. Oppenheim, “Advanced Topics in Signal Processing,” 
Prentice Hall: Englewood Cliffs, NJ, 1988. 

[9] C. S. Burrus, R. A. Gopinath, H Guo, Introduction to Wavelets and Wavelet 
Transforms A Primer, Prentice Hall: Upper Saddle River, NJ, 1998. 

[10] G. H. Watson, “Introduction to Wavelet Analysis,” NATO RTO Lecture Series 
216, Application of Mathematical Signal Processing Techniques to Mission 
Systems, 1999. 

[11] G. H. Watson, “The Detection of Unusual Events in Cluttered Natural 
Backgrounds,” NATO RTO Lecture Series 216, Application of Mathematical 
Signal Processing Techniques to Mission Systems, 1999. 

[12] G. Strang, Linear Algebra and Its Applications, Harcourt Brace Jovanovich 
College Publishers: New York, NY, 1988. 

[13] K. Ogata, State Space Analysis of Control Systems, Prentice Hall: Englewood 
Cliffs, NJ, 1967. 



171 



[14] J. Burl, Linear Optimal Control, Addison Wesley: Berkely, CA, 1999. 

[15] M. Misti, Y. Misti, and G. Oppenheim, Wavelet Toolbox for Use with MATLAB, 
The Math Works Inc.: Natick, MA, 1997. 

[16] J. Proakis, D. Manolakis, Introduction to Digital Signal Processing, Macmillan: 
New York, NY, 1988. 



172 



INITIAL DISTRIBUTION LIST 



1. Defense Technical Information Center 2 

8725 John J. Kingman Rd„ STE 0944 

Ft. Bel voir, VA 22060-6218 

2. Dudley Knox Library 2 

Naval Postgraduate School 

41 1 Dyer Rd. 

Monterey, CA 93943-5101 

3. Mr. John Cannon 1 

Advanced Deployable System, PMW 183-6 

Space and Naval Warfare Systems Command (SPAWAR) 

4301 Pacific Highway (OTI) 

San Diego, CA 92110-3127 

4. Chairman, Code EC 1 

Department of Electrical and Computer Engineering 

Naval Postgraduate School 
Monterey. CA 93943-5121 

4. Prof. Roberto Cristi, Code EC/Cx : 2 

Department of Electrical and Computer Engineering 

Naval Postgraduate School 
Monterey, CA 93943-5121 

5. Prof. Monique P. Fargues, Code EC/Fa 1 

Department of Electrical and Computer Engineering 

Naval Postgraduate School 
Monterey, CA 93943-5121 

6. Mr John A. Thornton 1 

Chief Engineer 

Advanced Deployable System, PMW 183-3 

Space and Naval Warfare Systems Command (SPAWAR) 

4301 Pacific Highway (OTI) 

San Diego, CA 92110-3127 

7. Robert E. Stevens 1 

5601 Bideford Ct 

Bowie. MD 20715 



173 



8. LT John D. Stevens 3 

1 16 E. Berger St 
Emmaus, PA 18049 



174 



60 “THU 

6/02 22527 -200 >* 6 



